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1  Summary 

The  extreme  impact  loading  consider  in  this  work  consists  of  high  amplitude 
loads  with  short  durations  which  can  result  in  high  frequency  content.  These 
conditions  can  result  in  wave  propagation  with  significant  nonlinear  character¬ 
istics.  With  the  use  of  conventional,  linear  high  fidelity  methods,  the  effects 
of  the  nonlinearity  can  be  overlooked  and  its  inffuence  on  the  structure  will 
not  be  known.  The  focus  of  this  work  was  the  development  of  high  fidelity 
modeling  and  analysis  methods  for  studying  this  nonlinear  behavior.  The 
Alternating  Frequency-Time  Finite  Element  Method  (AFT-FEM)  was  devel¬ 
oped  from  the  spectral  finite  element  method  (SFEM)  in  order  to  expand  its 
high  fidelity  performance  to  study  nonlinear  wave  propagation.  Later  efforts 
produced  the  Alternating  Wavelet-Time  Finite  Element  Method  (AWT-FEM) 
by  using  a  wavelet  basis.  A  wavelet  basis  is  better  conditioned  for  study¬ 
ing  localized  and  transient  signal  characteristics  such  as  those  encountered  in 
impact  loading  and,  as  a  result,  AWT-FEM  provides  further  improvements 
beyond  AFT-FEM.  A  practical  force  identification  scheme  was  also  developed 
by  using  SFEM  to  determine  both  temporal  and  spatial  information  about  an 
impact  load  by  using  only  a  small  number  of  sensors.  The  force  identification 
scheme  was  first  applied  to  linear  systems  and  then,  by  using  the  concepts 
developed  for  AWT-FEM,  it  was  extended  to  be  able  to  perform  force  identifi¬ 
cation  for  nonlinear  systems.  These  numerical  simulation  techniques  and  force 
identification  methods  provide  the  means  to  gain  new  insight  into  the  effects 
of  extreme  impact  loading  on  structures.  This  improved  understanding  can 
them  be  used  to  design  new  systems  and  structures  and  to  better  understand 
how  they  will  perform  under  severe  conditions. 

Key  words:  Impact  Response,  Nonlinear  Wave  Propagation,  Spectral  Finite 
Element  Method,  Wavelets,  Force  Identification 


2  Introduction 


In  the  majority  of  traditional  scientific  and  technological  applications,  wave 
propagation  in  solid  structures  is  bounded  by  the  elastic  limit  and  a  linear 
approximation  is  often  adopted.  However,  as  the  exploration  is  entering  into 
more  advanced  helds  like  new  materials,  extreme  working  environments,  and 
high  hdelity  analysis,  nonlinearity  can  be  introduced  in  various  ways  and  can 
become  a  crucial  factor  that  must  not  be  neglected.  In  the  past  decade,  fa¬ 
cilitated  by  the  development  of  numerical  techniques  and  high  performance 
computers,  signihcant  progress  has  been  made  in  the  study  of  nonlinear  elas¬ 
tic  wave  propagation  for  applications  including  biomedical  analysis  [1],  non¬ 
destructive  evaluation  [2],  and  seismic  motion  analysis  [3]. 

In  many  military  and  civil  applications,  nonlinear  mechanical  wave  prop¬ 
agation  occurs  through  structural  components  that  can  be  approximated  as 
rods,  beams,  or  plates.  For  example,  in  the  force  identihcation  of  aircraft 
through  impact  wave  analysis,  the  impact  load  often  has  a  high  amplitude 
with  a  short  duration,  where  high  frequency  content  in  the  response  will  sig- 
nihcantly  strengthen  the  influence  of  the  nonlinear  properties.  The  wings  of 
the  aircraft  can  be  approximated  as  a  beam  or  plate  structure  with  a  geometric 
nonlinearity.  Another  example  is  in  the  monitoring  and  control  of  drill  strings 
in  the  oilheld.  Drill  strings  can  be  modeled  as  long  geometrically  nonlinear 
beams  with  different  cross  sections  and  material  properties  [4].  Elastic  waves 
caused  by  impacts  with  the  wellbore  and  stress  waves  in  the  mud  flow  can 
be  used  to  control  the  stability  [5].  A  third  example  is  in  the  dynamics  of 
high-speed  and  high-precision  parallel  robots  [6].  Each  link  of  the  robot  can 
be  modeled  as  a  beam  undergoing  large  deformation.  For  these  systems  and 
many  other,  exact  solutions  of  linear  wave  equations  for  rod,  beam,  and  plate 
structures  are  known  [7]  and  analytical  relationships  in  the  nonlinear  wave 
equation  are  also  well  established  [8].  Incorporating  those  results  with  a  gen¬ 
eral  numerical  framework  will  greatly  improve  the  computational  performance 
and  provide  signihcantly  improved  accuracy. 

For  these  systems  and  many  other,  exact  solutions  of  linear  wave  equations 
for  rod,  beam,  and  plate  structures  are  known  [7]  and  analytical  relationships 
in  the  nonlinear  wave  equation  are  also  well  established  [8].  Incorporating 
those  results  with  a  general  numerical  framework  will  greatly  improve  the 
computational  performance  and  provide  signihcantly  improved  accuracy. 


2.1  Modeling  Nonlinear  Wave  Propagation 

Based  on  the  specific  nature  of  nonlinear  elastic  wave  propagation,  a  good 
numerical  technique  should  fulfill  three  requirements.  (1)  Provide  the  means 
to  predict  high  fidelity  responses  -  Nonlinear  phenomenon  can  be  accurately 
captured  in  the  response  and  available  for  further  analysis.  (2)  Provide  high 
computational  efficiency  -  The  method  should  take  advantage  of  current  com¬ 
putation  technology  to  ensure  studying  larger  models  is  not  computational 
prohibitive.  (3)  Be  flexible  -  The  method  should  be  able  to  be  applied  to 
various  structures  for  a  range  of  working  constraints. 

A  comparison  of  existing  methods  is  shown  in  Table.  1.  There  are  five 
examination  criteria:  (1)  high  fidelity  performance,  (2)  applicability  to  study 
linear  behavior,  (3)  applicability  to  study  nonlinear  behavior,  (4)  ability  to 
represent  physically  realistic  boundary  constraints,  and  (5)  parallel  computing 
compatibility.  The  check  mark  (/)  indicates  that  the  criterion  is  satisfied,  the 
dash  mark  (— )  indicates  the  criterion  is  partially  satisfied,  and  the  X  mark  (x) 
indicates  that  the  criterion  is  not  satisfied.  The  AFT-FEM  and  AWT-FEM 
developed  in  this  work  are  included  for  comparison. 


Table  1:  Comparison  of  different  methods. 


Methods 

High  Fidelity 

Linear 

Nonlinear 

Boundary 

Parallel 

TFEM  [9] 

X 

/ 

/ 

/ 

— 

MTFEM  [10] 

/ 

/ 

/ 

/ 

— 

SFEM  [11] 

/ 

/ 

X 

— 

/ 

WSFEM  [12] 

/ 

/ 

X 

/ 

/ 

AFT-FEM 

/ 

/ 

/ 

— 

/ 

AWT-FEM 

/ 

/ 

/ 

/ 

/ 

The  finite  element  method  (FEM)  [9]  has  been  widely  applied  to  study 
structural  dynamics  and  has  developed  into  a  class  of  numerical  methods  with 
various  time  integration  and  space  discretization  techniques.  In  Table  1,  the 
time-domain  FEM  (TFEM)  refers  to  a  standardized  procedure  composed  of  a 
serial  time-integration  technique  like  the  Newmark  method,  an  iteration  ap¬ 
proach  like  the  Newton-Raphson  method  to  address  nonlinearity,  and  a  stan¬ 
dard  finite  element  space  discretization  to  approximate  the  geometry.  Other 
modified  versions  of  FEM  with  various  time-integration  techniques  like  space- 
time  finite  element  method  [13]  is  not  considered  here.  TFEM  can  be  used 
to  solve  both  linear  and  nonlinear  problems  with  physically  realistic  boundary 
conditions  and  complex  geometries.  However,  for  transient  wave  propagation. 


especially  those  caused  by  impact  loading  with  high  frequency  content,  re¬ 
sponses  predicted  with  this  method  can  exhibit  spurious  oscillations  which  re¬ 
sults  from  erroneously  introduced  numerical  dispersion  and  dissipation  [10,  14], 
This  error  will  accumulate  as  the  wave  travels  along  a  structure  and  signif¬ 
icantly  distort  the  response  unless  an  extremely  hne  meshing  is  adopted  to 
capture  the  sharp  gradient  change  in  the  wave  shape.  Also,  due  to  the  serial 
nature  of  the  Newmark  method,  TFEM  cannot  be  directly  incorporated  into 
a  parallel  computing  framework.  Speed-up  procedures  are  needed  to  make  it 
compatible  with  a  parallel  computational  framework. 

In  order  to  solve  the  spurious  oscillation  error  and  obtain  a  high  hdelity 
response,  many  modihed  time-domain  hnite  element  methods  (MTFEM)  have 
been  developed  by  using  more  advanced  interpolation  functions  [10,  15].  The 
spectral  element  method  (SEM)  Erst  proposed  by  Patera  [16]  in  fluid  dynamics 
has  been  adapted  to  study  ID  linear  wave  propagation  in  solid  structures  [15]. 
In  SEM,  high  order  polynomials  have  been  chosen  as  local  shape  functions  and 
a  high  degree  of  accuracy  in  the  response  can  be  observed.  An  application  of 
SEM  to  study  ID  nonlinear  wave  propagation  in  a  rod  was  reported  in  [17]. 
Even  though  theoretically  SEM  can  be  adapted  to  solve  nonlinear  problems, 
other  related  works  are  surprisingly  rare  in  the  literature.  Ham  and  Bathe 
proposed  another  modihed  method  using  low-order  polynomials  enriched  with 
harmonic  functions  [10].  A  specihc  scheme  has  also  been  designed  to  overcome 
ill-conditioning.  This  method  can  be  used  to  solve  both  linear  and  materially 
nonlinear  systems,  but  may  encounter  challenges  for  systems  with  geometric 
nonlinearities.  Extra  efforts  on  reformation  are  needed  to  make  these  methods 
compatible  for  the  parallel  computing  framework. 

Another  class  of  methods  transforms  the  model  and  solution  procedure 
into  a  spectral  domain.  A  typical  example  is  the  fast  Fourier  transform  (FFT)- 
based  spectral  finite  element  method  (SFEM)  proposed  by  Doyle  [7,  11].  In  the 
literature,  there  is  some  confusion  over  the  terms  SFEM  and  SEM  [15].  Here, 
SFEM  specihcally  refers  to  the  Fourier-based  method  developed  by  Doyle.  For 
linear  problems,  SFEM  uses  exact  wave  solution  as  frequency-dependent  shape 
functions  and  only  one  element  is  required  to  obtain  highly  accurate  responses 
for  a  given  position  in  a  structural  component  like  a  rod  or  a  beam.  For  a 
structure  like  a  frame  or  a  truss,  an  element-based  procedure  can  be  adopted 
at  each  independent  frequency.  This  makes  the  method  inherently  suitable 
for  parallel  computing.  With  the  incorporation  of  the  exact  analytical  wave 
solutions,  the  capacity  for  parallel  computing  and  the  fast  computational  per¬ 
formance  of  FFT,  the  SFEM  can  be  utilized  as  a  real-time  method  to  obtain 
high  hdelity  responses  to  impacts  with  high  frequency  content.  Various  appli¬ 
cations  on  linear  ID  and  2D  problems  can  be  found  in  Gopalakrishnan  and 


Lee’s  works  [18,  19].  However,  this  method  has  two  major  drawbacks.  First, 
the  periodic  nature  of  the  basis  functions  used  for  the  Fourier  transform  can 
result  in  wrap-around  issues  in  which  the  signal  outside  the  time-window  will 
become  wrapped  around  to  the  beginning  of  the  signal  sequence.  Semi-infinite 
elements  (non-reflecting  boundary  condition)  are  often  added  at  the  ends  of  a 
structure  to  leak  energy  out  of  the  system.  For  other  types  of  boundary  con¬ 
ditions,  additional  damping  combined  with  a  longer  time-window  is  required. 
This  greatly  deteriorates  the  computational  performance  of  SFEM  when  it 
is  applied  to  study  wave  propagation  in  structures  with  physically  realistic 
boundary  constraints.  Second,  nonlinear  terms  often  result  in  convolution 
terms  in  the  frequency-domain,  which  are  computationally  prohibitive  to  cal¬ 
culate  with  iterative  procedures.  This  makes  SFEM  not  directly  applicable  for 
nonlinear  systems. 

In  order  to  avoid  the  wrap-around  issue  of  SFEM,  Mitra  and  Gopalakrish- 
nan  developed  a  wavelet-based  spectral  hnite  element  method  (WSFEM)  [12]. 
Instead  of  using  FFT,  a  discrete  wavelet  transform  (DWT)  on  the  sampling 
scale  was  adopted  to  transform  the  linear  wave  equation  into  a  wavelet-domain 
where  each  wavelet  point  represents  a  shifting  along  time  axis.  A  wavelet  ex¬ 
trapolation  technique  was  also  incorporated  to  represent  non-periodic  bound¬ 
ary  conditions  [20].  WSFEM  successfully  circumvented  the  wrap-around  issue 
and  can  also  be  incorporated  into  a  parallel  computing  framework  after  being 
uncoupled  through  an  eigenvalue  analysis.  Applications  to  linear  wave  propa¬ 
gation  in  composite  and  nano-composite  structures  were  reported  in  [21].  An 
adaptation  of  the  wavelet  transformation  technique  in  WSFEM  with  SEM  spa¬ 
tial  interpolation  was  also  applied  to  2D  and  3D  linear  elastic  waves  [22,  23]. 
Similar  to  SFEM,  WSFEM  cannot  be  directly  applied  to  study  nonlinear  sys¬ 
tems. 

An  intuitive  solution  to  deal  with  the  nonlinear  problem  in  frequency- 
domain  methods  is  to  transform  the  calculation  of  the  nonlinear  terms  back  to 
the  time-domain.  An  alternating  frequency-time  iterative  framework  is  needed 
to  transform  the  result  of  the  nonlinear  calculation  into  the  frequency-domain, 
solve  for  the  linear  solution  by  using  SFEM,  and  then  return  to  the  time- 
domain  to  calculate  new  values  of  the  nonlinear  terms  for  the  next  iteration. 
This  framework  was  hrst  proposed  by  Cameron  [24].  A  combination  of  it  with 
the  SFEM  was  performed  by  Lee  [25]  to  study  ID  nonlinear  blood  flow  in 
human  arteries. 


2.2  Force  Identification 


Indirect  methods  for  impact  force  identification  have  attracted  researchers 
dne  to  the  nonlinearity  of  the  impact  problem  and  complexity  of  impact  inci¬ 
dents  [26,  27,  28,  29].  Nnmerous  techniques  have  been  developed  which  uses 
inverse  methods  for  impact  force  identihcation.  Some  common  techniques  for 
force  identihcation  are  the  deconvolution  method  [26,  27],  state  variable  for¬ 
mulation  [30],  the  sum  of  weighted  accelerations  [31],  and  the  spectral  hnite 
element  method  [12],  The  most  popular  method  is  the  deconvolution  tech¬ 
nique  which  uses  the  assumption  of  a  linear  response  in  order  to  allow  for  the 
application  of  the  convolution  integral  to  identify  the  force  information.  For 
linear  systems,  the  convolution  integral  of  a  system’s  impulse  response  and  the 
applied  force  results  in  the  response  of  the  system.  This  technique  has  been 
applied  in  both  the  time-domain  [26]  and  the  frequency- domain  [32]. 

In  a  study  by  Doyle,  a  time-domain  deconvolution  technique  was  success¬ 
fully  developed  in  order  to  experimentally  obtain  dynamic  contact  laws  [26]. 
Response  behavior  was  monitored  by  using  strain  gages  affixed  to  a  beam 
structure  and  the  impulsive  load  was  applied  by  using  a  pendulous  ball.  In 
other  work,  Chang  and  Sun  calculated  the  applied  impact  force  by  using  an 
experimental  Green’s  function  and  time-domain  signal  deconvolution  [33] .  The 
reconstructed  force  was  identified  independent  of  the  location  of  the  sensors 
on  a  composite  beam  structure.  In  both  cases,  the  position  of  the  impact  was 
assumed  to  be  known. 

The  deconvolution  technique  is  also  easily  implemented  in  the  frequency- 
domain  for  conducting  force  identification  due  to  the  ease  of  performing  the 
deconvolution  calculation  since  the  convolution  integral  in  the  time-domain 
corresponds  to  a  simple  multiplication  operation  in  the  frequency-domain.  In 
other  work  of  Doyle,  a  frequency-domain  deconvolution  method  was  devel¬ 
oped  and  used  to  experimentally  determine  contact  laws  [32].  The  impact 
force  was  estimated  by  using  analytical  relationship  between  contact  force  and 
the  strain  response  of  structure.  In  order  to  improve  the  quality  of  the  recon¬ 
structed  force,  signal  processing  techniques  were  applied  on  the  strain  response 
to  avoiding  the  wrap-around  problem  [34,  35,  36]. 

Similar  to  the  force  identification  methods,  indirect  methods  of  identifica¬ 
tion  of  the  impact  location  have  been  studied  extensively  due  to  limitations  in 
the  direct  identification  of  the  impact  location  [26,  32,  37,  38,  39].  Methods 
for  the  localization  of  impacts  can  be  classified  into  two  main  categories.  The 
first  type  of  method  uses  wave  propagation  speed  and  arrival  time  to  determine 
impact  location  [40,  41,  42,  43,  44].  Methods  in  the  first  category  take  advan¬ 
tage  of  the  wave  speed  and  the  recorded  time  for  the  traveling  waves  to  reach 


the  sensors  in  order  to  calculate  the  position  of  the  impact  force.  Gaul  and 
Hurlebaus  developed  a  method  based  on  the  arrival  times  of  waves  at  various 
frequencies  [40] .  The  location  of  the  impact  was  identihed  by  solving  a  nonlin¬ 
ear  equation  extracted  from  the  speed  of  the  traveling  waves.  The  second  type 
of  localization  of  impact  force  is  performed  by  minimizing  an  error  function 
which  is  dehned  based  on  the  impact  force  and  the  system  response  [45,  46,  47]. 
The  system  response  can  be  obtained  from  a  system  model  or  transfer  function 
which  is  derived  analytically,  numerically,  or  experimentally  [48]. 

Additionally,  it  is  possible  for  the  two  methods  to  be  used  simultaneously  to 
accurately  locate  the  impact  force.  In  the  recent  work  of  Liang  et  al.,  elements 
of  the  two  methods  were  combined  to  create  a  distributed  coordination  algo¬ 
rithm  for  locating  the  impact  force  with  greater  computational  efficiency  [49] . 
However,  further  studies  are  required  in  order  to  fully  explore  the  effectiveness 
of  this  approach. 

2.3  Overview 

In  the  section  3,  the  key  elements  of  the  two  techniques  for  (i)  modeling  non¬ 
linear  wave  propagation  and  (ii)  performing  force  identihcation  are  presented. 
The  performance  of  these  methods  is  demonstrated  in  section  4.  The  main  con¬ 
clusions  of  this  work  and  the  direction  of  future  work  are  presented  in  section  5. 
The  material  presented  in  this  report  constitutes  the  main  advances  which 
resulted  from  this  work.  For  full  details,  refer  to  the  publication  produced 
through  this  work  [50,  51,  52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  64,  65]. 

3  Methodology 

In  this  section,  the  methods  developed  in  this  work  are  presented.  First,  the 
Alternating  Wavelet-Time  Finite  Element  Method  (AWT-FEM)  is  presented. 
Then,  the  force  identihcation  and  force  localization  methods  are  presented. 

3.1  Modeling  Nonlinear  Wave  Propagation 

A  specihc  wavelet  transform  is  described  in  the  hrst  subsection.  The  spectral 
hnite  element  method  formulation  in  the  wavelet- domain  is  discussed  in  the 
second  section.  The  procedure  of  AWT-FEM  is  described  in  the  third  subsec¬ 
tion.  Its  applications  to  an  elementary  rod  model  with  a  material  nonlinearity, 
an  Euler-Bernoulli  beam  model  with  a  geometric  nonlinearity,  and  a  conven¬ 
tional  plate  model  with  a  weak  geometric  nonlinearity  are  derived  in  the  last 


subsection. 


3.1.1  Spectrally-Uncoupled  Wavelet  Transform 

The  general  governing  differential  wave  equation  of  a  nonlinear  structure  is 
given  as  [25] 


Tu(x,  t)  +  Alu(x,  t)  +  AAu(x,  t)  +  f  (u,  X,  t)  =  p(x,  t),  (1) 

where  C  is  the  linear  differential  operator  with  respect  to  the  spatial  coordinate 
vector  X,  the  over-dot  represents  a  derivative  with  respect  to  time  t,  M  is  the 
inertial  operator,  Af  is  the  damping  operator,  u(x,  t)  is  the  displacement  held 
vector,  p(x,  t)  is  the  external  load  vector,  and  f(x,  t)  is  the  nonlinear  term 
vector. 

A  single-scale  wavelet  Galerkin’s  method  using  a  Daubechies  wavelet  is 
applied  to  transform  Eqn.  (1)  into  a  set  of  spectrally- uncoup  led  partial  differ¬ 
ential  eqnations  (PDEs).  Details  of  this  method  and  its  application  to  a  ID 
linear  wave  equation  can  be  fonnd  in  [12,  20,  66].  Here,  the  application  of  this 
method  to  the  general  nonlinear  wave  equation  in  Eqn.  (1)  is  briehy  described 
in  the  following  three  steps. 

Step  1: 

A  new  variable  r  =  t/At  is  introduced.  In  the  following  derivation,  r  is 
regarded  as  a  continuous  variable  since  integration  with  respect  to  r  will  be 
performed. 

u(x,  t)  is  represented  in  an  expansion  form  by  translating  the  wavelet  scal¬ 
ing  function  <yc(r). 

u(x,t)  =  u(x,t)  = '^Uk(x)(f(T  -  k),  kez,  (2) 

k 

where  Ufc(x)  (referred  as  hereafter)  is  the  approximation  coefficient  of  u(x,  t) 
at  an  associated  point  k,  Z  is  the  set  of  integers..  By  utilizing  the  orthogonality 
of  the  scaling  function,  the  approximation  coefficient  is  obtained  as 

Ufc  =  At  J  u{x,T)ip{T  —  k)dT.  (3) 

The  scaling  function  (p{t)  is  dehned  on  the  hnest  scale  (the  sampling  scale). 
Hence,  the  integer  k  presents  a  shifting  along  the  sampled  time  axis  for  k  sam¬ 
ple  points.  It  is  comparable  to  the  frequency  component  at  a  given  frequency 
in  the  Fourier  expansion,  in  the  sense  that  both  the  expansion  in  Eqn.  (2)  and 


Fourier  expansion  expand  the  original  function  by  a  class  of  basis  function 
indexed  by  a  variable.  The  difference  here  is  the  index  k  indicates  a  shifting 
in  time  but  the  corresponding  variable  in  Fourier  expansion  represents  a  re- 
hnement  of  frequency.  The  physical  meaning  of  k  indicates  that  the  expansion 
in  Eqn.  (2)  is  in  nature  a  shifting  tool  for  a  time  integration  technique.  The 
index  k  is  referred  as  “wavelet  point”  to  distinguish  the  time  series  index  r 
before  the  expansion.  Similarly,  the  nonlinear  term  f(x,  f)  and  the  loading 
term  p(x,  f)  are  expanded  by  ^k  and  as  follows. 


f(x,f)  =  ^ffc(p(r  - /c),  fc  e  Z, 

k 

(4) 

p(x,  t)  =  ^  Pk^{T  -k),  k  ez. 

(5) 

k 


The  above  expansions  are  substituted  into  Eqn.  (1).  By  multiplying  the 
resulting  equation  by  the  scaling  function  (p{T—j)  and  integrating  it  over  time, 
the  following  set  of  equations  are  obtained. 

j  -  k)^{T  -  j)dT+ 

J —  -  k)ip{T  -  j)dT  +  ij  =  pj,  j  =  0,l,...,n- 1,  (6) 

where  f(x, f),  p(x, t),  and  the  u(x, t)  after  the  linear  differential  operator  C 
are  transformed  into  their  corresponding  approximation  coefficients  fj ,  pj ,  and 
Uj  due  to  the  orthogonality  of  the  wavelet  scaling  function.  For  the  time 
derivatives  of  u(x,  r),  the  two  integrals  associated  with  Ai  and  Af  in  Eqn.  (6) 
are  evaluated  in  the  form  of  connection  coefficients  for  compactly  supported 
wavelets  [66],  as  dehned  in  Eqn.  (7)  and  (8).  The  index  j  represents  a  given 
wavelet  point.  Similar  to  the  dehnition  of  k,  it  originally  indicates  a  shifting 
along  the  time  axis. 

n]_k  =  j  0(r  -  k)^{T  -  j)dT,  (7) 

=  j  k)ip(T  -  j)dT.  (8) 

For  a  multidimensional  problem,  u  =  [u^v^wY  and  Uj  =  [uj,Vj,WjY ,  a 
complete  matrix  form  of  Eqn.  (6)  with  respect  to  the  three  dimensions  and  the 
wavelet  point  j  is  cumbersome.  The  couplings  between  different  dimensions 
{x,y,z)  are  assumed  to  be  represented  in  the  nonlinear  term  fj.  Since  fj  will 


be  treated  as  an  equivalent  nodal  force  term  without  explicitly  evaluating  it, 
and  only  its  nodal  values  will  be  used  in  the  numerical  procedure  to  solve  the 
equation,  Eqn.  (6)  can  be  regarded  as  decoupled  with  respect  to  dimensions. 
The  coupling  in  Eqn.  (6)  is  with  respect  to  the  index  j  and  is  captured  in  the 
matrix  f2.  Hence,  Eqn.  (6)  describes  a  set  of  spectrally  coupled  PDEs  with 
respect  to  j,  and  has  to  be  solved  for  each  wavelet  point  (translation  along 
time  axis)  j. 

Step  2: 

Next,  a  wavelet  extrapolation  technique  [20]  is  adopted  to  address  the  time 
boundaries.  Specihc  treatments  of  periodic  and  non-periodic  time  boundaries 
for  linear  wave  equations  can  be  found  in  [12].  It  uses  a  polynomial  of  order 
A^/2  — 1  to  extrapolate  the  values  at  the  boundaries  from  the  known  coefficients 
near  the  boundaries.  The  connection  coefficient  matrix  Q,  becomes  a  n  x  n 
matrix  A,  which  is  a  version  with  non-periodic  time  boundaries.  The  result 
of  this  process  is  arranged  into  three  decoupled  matrix  forms  with  respect  to 
the  three  dimensions,  each  of  which  will  be  further  decoupled  with  respect  to 
the  spectral  index  j  in  the  next  step. 

Step  3: 

A  standard  eigenvalue  analysis  is  performed  to  uncouple  the  matrix  form 
of  Eqn.  (6)  with  respect  to  j.  The  eigenvalue  and  the  nxn  eigenvector  matrix 
of  the  nxn  connection  coefficient  matrix  A  are  Aj,  j  =  0, 1,  •  •  • ,  n  —  1  and  <h, 
respectively.  After  each  dimension  is  spectrally  uncoupled,  the  result  can  be 
organized  as 


C\ij  +  [MX]  MXj)  Uj  ij  =  Pj,  j  =  0,1, . . .  ,n  -  1,  (9) 

where  the  eigenvalue  Xj  can  be  treated  as  a  pseudo-frequency  and  is  related  to 
the  physical  angular  frequency  ujj  in  the  Fourier  transform  through  Xj  =  iuj, 
and 
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where  the  bracket  represents  a  vector  formation  with  respect  to  j. 

These  three  steps  dehne  a  special  single-scale  wavelet  transform  to  reduce 
a  temporal  problem  into  a  set  of  spectrally-uncoupled  PDEs.  This  procedure 
involves  several  wavelet-based  numerical  techniques  [20,  66],  an  eigenvalue 
analysis,  and  an  inverse  calculation  of  matrices.  However,  it  is  worth  not¬ 
ing  that  this  procedure  is  independent  of  the  system  to  be  studied  and  must 
only  be  performed  once  for  a  given  wavelet  scaling  function  ip{t)  and  length 
of  time  window  n.  The  resulting  eigenvalues  Aj,  eigenvector  matrix  <h,  and 
the  inverse  of  <h  can  then  be  calculated  and  stored  for  future  use.  The  scaling 
function  is  dehned  at  the  sampling  rate  of  the  signal,  which  is  also  the  finest 
scale  achievable  for  a  given  set  of  data.  In  the  above  wavelet  transform  proce¬ 
dure,  the  approximation  coefficient  of  a  given  signal  is  hrst  obtained  by  using 
a  numerical  applying  of  Eqn.  (3).  The  vector  of  approximation  coefficients, 
which  has  the  same  length  as  the  original  signal,  is  multiplied  by  the  inverse 
of  the  eigenvector  matrix  $  in  order  to  obtain  the  corresponding  values  in  the 
dehned  wavelet-domain.  The  inverse  transform  is  then  performed  by  multipli¬ 
cation  by  the  eigenvector  matrix  <h  to  obtain  the  approximation  coefhcients  in 
the  time-domain.  The  signal  is  recovered  by  hnally  using  Eqn.  (2). 


3.1.2  Spectral  Finite  Element  Method  Formulation 

Various  methods  are  available  to  solve  the  linear  version  of  Eqn.  (9).  In  [22,  23], 
the  SEM  is  chosen  to  solve  general  linear  2D  and  3D  problems  in  the  wavelet- 
domain.  In  this  study,  the  research  focus  is  to  develop  a  numerical  technique 
for  structures  that  can  be  modeled  with  rod,  beam,  and  plate  elements.  Uti¬ 
lizing  the  analytical  wave  solutions  of  those  models  may  greatly  improve  the 
computational  efficiency  and  alleviate  potential  convergence  issues  with  the 
iterative  approach  that  will  be  constructed  to  address  nonlinearities. 

The  linear  homogeneous  wave  equation  is  obtained  by  removing  the  non- 


linear  term  and  external  force  term  in  Eqn.  (9). 

Cuj  +  (MXj  +  A/’Aj)  Uj  =  0,  j  =  0, 1, . . . ,  n  —  1.  (11) 

For  ID  problems,  x  =  x  and  Eqn.  (9)  is  a  set  of  nonlinear  ODEs.  Exact 
solutions  to  Eqn.  (11)  for  rod  and  beam  structures  are  well-established  [7]. 
In  linear  SEEM,  those  exact  solutions  are  directly  employed  to  solve  for  the 
response  at  a  given  location  x  and  only  one  element  is  needed  for  the  formula¬ 
tion.  In  this  nonlinear  derivation,  the  exact  linear  solutions  are  adopted  as  the 
shape  functions  in  a  standard  finite  element  framework  in  the  wavelet-domain. 
The  spectral  displacement  vector  iij  is  expressed  as 

Uj  =N{x-,Xj)qj,  (12) 

where  N(a:;  Xj)  (referred  as  N  hereafter)  is  the  spectral  shape  functions  based 
on  the  exact  linear  wave  solutions,  and  is  the  spectral  nodal  displacement 
vector  at  the  pseudo-frequency  Xj  or  wavelet  point  j.  By  following  a  stan¬ 
dard  Galerkin’s  method,  the  nonlinear  ODE  in  Eqn.  (9)  becomes  the  following 
matrix  form 

K(Aj)qj  =  qE  +  qN,  (13) 

where 


K(A,)  = 

J  (A(N,  N)  +  (M A^  +  AfXj)  N)  dx, 

qE  = 

/N-p,dx. 

qN  = 

(14) 

The  spectrally  formulated  dynamic  stiffness  matrix  is  represented  by  K(Aj), 
A(N,N)  is  the  function  resulting  from  integration  by  parts  [19],  qE  is  the 
equivalent  nodal  force  term  for  the  distributive  external  load,  and  qN  is  the 
equivalent  nodal  force  term  for  the  nonlinear  term.  It  is  worth  noting  that 
the  process  to  derive  Eqn.  (13)  may  involve  several  steps  of  iteration  by  parts, 
so  the  sign  on  the  right-hand  side  may  change  depends  on  the  explicit  form 
of  C.  The  expression  in  Eqn.  (13)  is  for  the  purpose  of  a  general  and  concise 
representation. 

For  a  2D  problem,  x  =  {x,  y}  and  Eqn.  (9)  is  a  set  of  nonlinear  PDFs. 
The  wavelet  transform  is  performed  again  on  Eqn.  (9)  with  respect  to  one  of 
the  spatial  coordinate  y.  After  that,  Eqn.  (9)  is  reduced  to  a  set  of  ODEs  and 
the  model  is  simplified  into  a  ID  problem.  A  similar  spectral  finite  element 


method  formulation  as  above  can  be  performed  to  obtain  a  matrix  form  as  in 
Eqn.  (13).  An  example  of  this  spatial  reduction  process  for  a  plate  model  is 
provided  at  the  end  of  this  section. 

3.1.3  General  Procedure  of  the  AWT-FEM 

Based  on  the  spectrally-uncoupled  wavelet  transform  and  the  spectral  hnite 
element  method  formulation  in  the  wavelet-domain,  an  alternating  wavelet¬ 
time  framework  is  constructed  to  iteratively  simulate  nonlinear  elastic  wave 
propagation  responses.  A  flowchart  of  the  AWT-FEM  is  presented  in  Fig.  1. 
The  initial  state  is  obtained  from  the  linear  model  by  setting  the  nonlinear 
term  f(u,  x,  t)  in  Eqn.  (1)  equal  to  zero.  The  f(x,  t)  is  written  as  f(u,  x,  f)  to 
indicate  that  the  value  of  the  nonlinear  force  term  is  directly  determined  by 
the  current  value  of  u.  At  each  iteration,  the  nodal  values  of  the  nonlinear 
term  f(u,  x,  f)  are  calculated  from  the  current  states  in  the  time-domain.  The 
result  and  the  nodal  values  of  external  force  p(x,  f)  are  transformed  into  the 
wavelet-domain  through  the  spectrally-uncoupled  wavelet  transform.  At  each 
pseudo-frequency  or  wavelet  point,  a  spectral  hnite  element  method  formula¬ 
tion  is  constructed  based  on  Eqn.  (13)  to  obtain  nodal  values  of  Uj.  Nodal 
values  of  the  derivatives  of  with  respect  to  x  are  easily  calculated  based  on 
Eqn.  (12).  Since  Eqn.  (9)  is  uncoupled  with  respect  to  j,  the  above  calculation 
is  performed  in  parallel  at  each  wavelet  point.  The  nodal  information  is  then 
transformed  back  to  the  time-domain  via  the  inverse  wavelet  transform.  New 
values  of  the  nonlinear  term  f(x,  f)  are  calculated  for  use  in  the  next  itera¬ 
tion.  The  process  is  continued  until  the  predehned  error  measure  dehned  in 
Eqn.  (15)  is  less  than  a  prescribed  tolerance. 

The  number  of  wavelet  point  is  determined  by  the  length  of  the  signal. 
The  higher  the  sampling  rate  and  the  longer  the  signal  length,  the  better 
accuracy  of  the  wavelet  transform  and  a  better  performance  of  the  AWT- 
FEM,  since  Eqn.  (2)  and  (3)  are  implemented  through  numerical  integration 
in  the  program.  In  most  of  the  applications,  the  signal  is  collected  through 
experiments  or  predehned  as  inputs.  Hence,  the  length  of  the  signal  or  the 
number  of  wavelet  point  is  taken  as  hxed  conditions  in  numerical  simulation 
in  this  study. 
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where  M  is  the  length  of  the  time  history,  xq  is  the  impacted  location,  and  the 
subscript  i  represents  the  iteration  number.  The  response  at  the  impacted  po- 


sition  has  both  the  initial  impulse  prohle  and  the  interaction  with  the  reflected 
wave  from  the  boundary  within  a  given  time-window.  It  can  reflect  most  of  the 
accumulated  nonlinear  information  in  the  wave  propagation  process.  Hence, 
the  impact  position  is  chosen  to  define  the  error  measure.  Other  point  might 
also  be  used  and  the  result  will  not  be  qualitatively  different. 
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Figure  1:  The  flowchart  for  the  AWT-FEM. 


3.2  Applications  to  Structural  Models 
(1)  Rod: 

A  rod  model  with  a  material  nonlinearity  is  considered  here.  A  general 
nonlinear  stress-strain  relationship  with  a  high  order  strain  term  has  been 
used  in  [8,  17,  67].  In  the  second  order  case,  the  nonlinear  constitutive  law 
becomes 

a  =  E{e  +  ae^),  (16) 

where  a  is  the  stress,  e  is  the  strain,  E  is  the  elastic  modulus,  and  a  is  the  non¬ 
linear  coefficient.  When  a  >  0,  the  quadratic  term  exhibits  a  hardening  effect 
on  the  constitutive  curve.  When  a  <  0,  a  softening  effect  can  be  observed. 

A  linear  strain-displacement  relationship  e  =  du/dx  is  assumed.  By  using 
Hamilton’s  principle,  the  governing  materially  nonlinear  homogeneous  wave 
equation  for  a  rod  is  obtained  as 


EAu”  —  pAu  +  E{x,  t)  =  0, 


(17) 


where 


F{x,t)  =  2aEAu'u'' ,  (18) 

the  cross  sectional  area  is  represented  by  A,  p  is  the  density,  u  is  the  axial 
displacement,  the  prime  represents  a  derivative  with  respect  to  the  spatial 
coordinate  x,  and  F[x,  t)  is  the  nonlinear  term. 

By  using  the  spectrally-uncoupled  wavelet  transform  previously  described, 
the  governing  equation  of  motion  in  the  wavelet-domain  is  obtained  as 

EAu'-  —  pAX'^Uj  +  Fj  =  0,  j  =  0, 1, . . . ,  n  —  1.  (19) 

The  expression  of  the  spectrally-dependent  shape  function  Nr  for  a  rod 
model  can  be  found  in  Eqn.  (3.36)  in  [19].  It  is  a  1  x  2  vector  with  com¬ 
ponents  as  functions  of  x  and  u.  The  angular  frequency  u  associated  with 
the  frequency-domain  implementation  in  the  references  must  be  replaced  by 
uj  =  —iXj  for  this  wavelet-domain  method,  which  is  indicated  by  comparing 
the  Eqn  (11)  in  the  wavelet-domain  with  Eqn  (3.4)  in  [19]  in  the  frequency- 
domain.  By  substituting  it  into  Eqn.  (14),  the  dynamic  stiffness  matrix  K  and 
the  equivalent  nodal  nonlinear  force  term  qn  become 

qn  =  y  ISi^Fjdx.  (20) 


(2)  Beam: 

A  beam  model  with  a  geometric  nonlinearity  is  considered  here.  The  strain 
is  assumed  to  be  small  and  the  rotation  is  expected  to  be  moderately  large.  The 
von  Kdrmdn  strain  [9]  is  dehned  to  represent  the  nonlinear  strain-displacement 
relationship.  The  constitutive  relationship  is  assumed  to  be  linear.  By  using 
Hamilton’s  principle,  the  governing  geometrically  nonlinear  homogeneous  wave 
equations  for  a  beam  are  obtained  as 


EAu”  —  pAu  +  Fji{x,  t)  =  0,  (21) 

EIw''"  +  pAiij  —  Fb{xA)  =  0,  (22) 

Fii{x^  t)  =  EAw'w",  (23) 

FB{x,t)  =  EAw"(u'+^{w'fy  (24) 


where 


the  axial  displacement  is  represented  by  m,  tn  is  the  transverse  displacement, 
I  is  the  second  moment  of  inertia,  Fji{x,t)  is  the  axial  nonlinear  force  term, 
and  Fsixjt)  is  the  transverse  nonlinear  force  term.  The  dehnitions  of  E  and 
A  are  same  as  those  for  the  rod  model. 

By  using  the  spectrally-uncoupled  wavelet  transform  developed  in  the  hrst 
subsection,  the  governing  equations  in  the  wavelet-domain  is  obtained  as 

EAv!'- -  pA\^-Uj  +  pRj^x)  =  0,  (25) 

Elw'j"  +  pAX^jWj  —  Fsjix)  =  0,  j  =  0, 1, . . . ,  n  —  1.  (26) 

The  expression  of  the  spectrally-dependent  axial  shape  functions  Nr  and 
transverse  shape  functions  Nr  can  be  found  in  Eqn.  (3.36)  and  Eqn.  (3.56)  in 
[19].  The  transverse  shape  function  Nr  is  a  1  x  4  vector  with  components  as 
functions  of  x  and  u.  As  was  done  for  rod  structures,  the  angular  frequency 
00  in  the  references  must  be  replaced  by  a;  =  —iXj-  By  organizing  Eqn.  (25) 
and  Eqn.  (26)  into  a  matrix  form  and  performing  the  spectral  hnite  element 
formulation,  the  dynamic  stiffness  matrix  K  becomes 
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where  Kmj  and  Ksij  are  components  of  Kr  and  Kr,  respectively. 
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(27) 


(28) 

(29) 


The  dynamic  stiffness  matrix  K  is  decoupled  between  axial  and  transverse 
motions.  The  equivalent  nodal  nonlinear  force  term  becomes 


/  NjupRjdx 
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(30) 


where  the  terms  Njn  and  Nbi  are  elements  of  Nr  and  Nr,  respectively.  The 


general  spectral  nodal  displacement  vector  is  dehned  as  = 
where  6  is  the  rotation  angle  in  the  wavelet-domain. 
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(3)  Plate: 

A  plate  model  with  a  geometric  nonlinearity  is  considered  here.  In  [9], 
the  classic  and  hrst-order  plate  theory  with  the  von  Kdrmdn  strain  is  chosen 
to  derive  the  governing  equation  of  motion.  The  results  are  three  coupled 
nonlinear  equations  in  the  three  spatial  dimensions  with  respect  to  two  in¬ 
plane  displacements  w,  n,  and  one  out  of  plane  displacement  w.  The  equation 
in  the  direction  perpendicular  to  the  plane  can  be  written  into  the  form  of 
von  Kdrmdn  equation  [68]  as 

DV^V‘^w{x,  y,  t)  +  phw{x,  y,  t)  =  F{x,  y,  t),  (31) 


where  =  d'^fdx'^  +  /dy"^  is  the  Laplace  operator,  h  is  the  thickness  of 

the  thin  plate,  w  is  the  transverse  displacement,  D  =  EK'/{12[1  —  is  the 
flexural  rigidity,  and  n  is  the  Poisson’s  ratio.  Other  parameters  are  dehned  as 
same  as  those  in  the  rod  model.  The  nonlinear  term  F{x,  y,  t)  is 


d‘^w 

F{x,  y,  t)  =  ^xx-^  +  N, 
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(32) 


where  Nyy^  and  N^y  are  resultant  membrane  forces  dehned  as 


rh/2 

rh/2 

rh/2 

1  ^yy 

1  (Jyydz^  d^xy 

1  ^xydz^ 

(33) 
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where  the  expressions  of  stresses  axxi  c^yy,  and  can  be  calculated  from  the 
strains  dehned  in  Eqn.  (6.2.5)  in  [9]. 

Equation  (32)  indicates  that  the  nonlinear  term  is  a  function  of  three  un¬ 
knowns  u,  V,  and  w,  which  correspond  to  displacement  along  x,  y,  and  2;  axis. 
Two  assumptions  are  made  here  to  further  simplify  the  model.  First,  it  is 
assumed  that  the  plate  is  only  subject  to  transverse  loading.  Second,  it  is  as¬ 
sumed  that  there  are  either  no  axial  constraints  on  the  boundaries  or  the  plate 
is  large  enough  so  that  the  propagating  wave  will  not  reach  the  boundaries 
within  the  time  window.  In  this  weak  geometric  nonlinear  case,  the  in-plane 
displacement  u  and  v  are  very  small  and  can  be  assumed  to  be  zero  in  Eqn.  (31) 
and  Eqn.  (32). 

By  using  the  spectrally- uncoup  led  wavelet  transform,  Eqn.  (31)  is  trans¬ 
formed  into  the  wavelet- domain  as 


DV^V'^Wj{x,y)  +  phX^-Wj{x,y)  =  Fj{x,y)  j  =  0, 1, . . . ,  n  -  1.  (34) 


The  spectrally-uncoupled  wavelet  transform  is  applied  again  with  respect 
to  one  of  the  spatial  dimension  y  and  Eqn.  (34)  is  reduced  to  a  ID  problem  as 


where  /3j  is  a  pseudo  spatial  frequency.  For  more  information  about  this  spa¬ 
tially  reduction  process,  refer  to  [21]  where  a  similar  example  for  a  linear  plate 
model  is  presented  (Chapter  7.3,  Eqn.  (7.34)  to  Eqn.  (7.37)). 

Equation  (35)  is  analogous  to  a  beam  equation.  By  setting  Fij  as  zero, 
Eqn.  (35)  becomes  a  linear  homogeneous  wave  equation.  The  exact  solution  is 
assumed  to  be  Wij  =  and  substituted  into  the  linear  homogeneous  wave 

equation.  Four  values  of  the  wavenumber  /c(/5j,  \j)  are  obtained  through  an 
eigenvalue  analysis.  The  spectrally-dependent  shape  function  Nb  is  dehned  in 
Eqn.  (3.56)  in  [19].  The  wavenumber  k  in  the  references  must  be  replaced  by 
the  new  calculated  values.  The  dynamic  stiffness  matrix  K  and  the  equivalent 
nodal  nonlinear  force  term  qN  are  derived  as 


(36) 


(37) 


where  explicit  expressions  of  integral  terms  can  be  referred  to  the  ones  in  the 
beam  model. 

3.3  Impact  Force  Identification 

In  the  impact  force  identihcation  method  presented  in  this  section,  the  dynamic 
stiffness  matrix  is  prepared  for  the  structure  and  multiplied  by  the  displace¬ 
ment  and  slope  information  in  order  to  calculate  the  identihed  force.  However, 
incomplete  deflection  and  slope  data  results  in  errors  in  the  calculated  force 
information.  Mechanical  wave  propagation  starts  immediately  after  the  im¬ 
pact  force  is  applied  and  it  moves  out  from  this  location  until  it  reaches  the 
ends  of  the  structure  and  is  reflected.  In  order  to  capture  the  complete  wave 
propagation  data  across  the  entire  beam  structure,  accelerometers  must  be 
mounted  along  the  entire  length  of  the  beam,  including  at  the  ends.  However, 
mounting  accelerometers  at  the  ends  of  the  structure  and  densely  along  the 
length  is  not  always  practical.  Therefore,  in  this  work  a  force  identihcation 


method  is  presented  that  utilizes  the  data  of  the  propagating  wave  recorded 
by  using  accelerometers  located  near  the  location  of  the  impact.  Although  the 
exact  location  of  the  impact  is  unknown,  it  is  assumed  that  the  general  area 
in  which  the  impact  occurs  is  known. 
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Figure  2:  Simulated  wave  propagation  (acceleration  data)  in  a  free-free  beam. 
The  boxed  area  corresponds  to  data  which  is  selected  for  analysis.  This  area 
corresponds  to  only  a  section  of  the  structure  and  does  not  include  reflections. 

The  simulation  acceleration  data  of  wave  propagation  through  the  beam 
structure  is  illustrated  in  Fig.  2.  The  hgure  illustrates  the  variation  of  acceler¬ 
ation  at  each  location  of  the  beam  through  time.  This  response  is  used  in  order 
to  identify  the  impact  force.  After  the  impact  occurs,  mechanical  waves  travel 
through  the  beam  structure  and  reflect  when  they  reach  the  ends.  The  reflec¬ 
tions  and  dispersion  in  beam  structures  increase  the  difficulty  of  performing 
successful  force  identihcation,  especially  if  it  is  not  possible  to  mount  sensors 
at  the  ends  of  the  structure.  In  this  force  identihcation  method,  the  impact 
force  is  calculated  by  using  a  STEM  model  and  a  subset  of  the  response  data. 
This  subset  of  the  displacement  and  slope  data  corresponds  to  a  reduction 
in  both  the  spatial  and  temporal  dimensions  around  the  location  where  the 
propagated  wave  originated,  as  illustrated  by  the  box  in  Fig.  2.  Note  that 
acceleration  data  is  presented  in  this  hgure.  This  subset  of  data  chosen  for  the 
identihcation  procedure  is  selected  so  that  it  does  not  contain  rehections.  The 
data  subset  is  used  with  a  beam  segment  model  to  calculate  the  desired  force 
information. 

A  howchart  of  the  force  identihcation  procedure  is  presented  in  Fig.  3.  The 


Figure  3:  Flowchart  illustrating  the  main  steps  of  the  force  identihcation 
scheme. 


acceleration  data  from  the  structure  is  transformed  into  the  frequency-domain 
by  using  the  fast  Fourier  transform  (FFT)  and  is  integrated  twice  to  obtain 
displacement.  The  slope  information  required  for  the  identihcation  process  is 
calculated  by  using  a  shape  function-based  method.  Zero-padding  is  applied 
to  the  response  data  in  order  to  address  the  FFT’s  assumption  of  periodicity. 
The  dynamic  stiffness  matrix  of  the  SFEM  model  for  a  segment  of  the  beam 
is  used  to  calculate  the  force  information  in  the  frequency-domain.  Detailed 
descriptions  of  each  step  are  presented  in  the  following  subsections. 

3.3.1  Time-Domain  Signal  Conditioning 

The  force  identihcation  procedure  uses  the  propagating  wave  data  immediately 
after  the  impact  to  identify  the  impulsive  force.  A  simple  algorithm  is  applied 
to  the  acceleration  data  in  order  to  identify  the  general  area  of  the  beam 
where  the  impact  occurred.  The  geometry  and  material  properties  of  the 
structure  are  used  to  approximate  the  length  of  time  after  the  impact  when 
rehections  are  expected  to  occur.  Exponential  windowing  is  applied  to  the 
complete  time  series  of  the  experimental  acceleration  data  in  order  to  reduce 
frequency-domain  leakage  in  the  response  [69].  Exponential  windowing  causes 
the  acceleration  to  decay  to  zero  at  the  end  of  the  full  data  set  and  ensures 
periodicity  in  the  response.  Since  only  a  small  portion  of  the  data  immediately 
after  the  impact  is  used,  this  windowing  has  a  negligible  effect  on  that  data 


and  subsequently  the  identified  force  information. 

Since  the  FFT  algorithm  assumes  periodicity,  it  can  result  in  wrap-around 
issues.  Zero  padding  is  later  applied  in  order  to  address  this  concern.  Zero 
padding  has  been  shown  to  significantly  improve  the  quality  of  force  informa¬ 
tion  identified  by  using  frequency- domain  methods  [32]. 

3.3.2  Deflection  Information 

Deflection  information  is  required  for  the  force  identification  procedure.  The 
deflections  of  the  nodes  are  calculated  from  the  acceleration  data.  The  inte¬ 
gration  is  performed  in  the  frequency-domain  in  order  to  improve  accuracy 
and  ensure  compatibility  with  the  frequency-domain  representation  of  the  sys¬ 
tem  [70].  The  integration  is  performed  by  dividing  the  frequency-domain  ac¬ 
celeration  data  by  the  imaginary  form  of  the  frequency.  For  the  zero  frequency, 
a  DC  value  is  used  for  the  integration  which  is  the  summation  of  the  acceler¬ 
ation  values  in  the  frequency-domain  [70].  However,  this  has  been  shown  to 
result  in  drift  in  the  time  history  of  displacement  data.  Fortunately,  the  small 
portion  of  the  displacement  data  right  after  the  impact  which  is  used  for  the 
force  identification  problem  is  negligibly  influenced  by  the  drift. 

3.3.3  Slope  Calculation 

The  acceleration  data  only  provides  information  regarding  the  one-dimensional 
translational  motion  at  discrete  points  along  the  structure  but  application  of 
the  force  identification  procedure  to  a  beam  structure  requires  both  deflection 
and  slope  information.  Therefore,  it  is  necessary  to  calculate  the  corresponding 
slope  response  at  each  location.  This  is  accomplished  by  taking  advantage 
of  the  kinematics  of  the  beam  structure.  In  the  slope  calculation  process, 
the  slope  information  of  the  beam  structure  is  obtained  by  using  the  shape 
function  of  the  SFEM  beam  model  of  the  structure.  This  technique  is  precise 
when  the  kinematics  of  the  structure  is  known.  In  order  to  calculate  the  slopes 
by  using  shape  functions,  the  beam  structure  is  meshed  with  four-node  beam 
elements  that  provide  continuity. 

The  slope  calculation  procedure  is  performed  in  the  frequency-domain  by 
using  deflection  data  at  each  node.  The  deflection  coefficients.  A,  B,  C,  and  D, 
of  each  four-node  beam  element  are  calculated  by  substituting  the  deflection 
values  at  the  four  nodes  for  the  specified  range  of  frequencies  into  the  deflec¬ 
tion  equation  and  solving  for  the  deflection  coefficients.  Slope  information  is 
obtained  by  using  the  analytical  derivative  of  the  shape  functions. 
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Figure  4:  Simulated  mechanical  wave  propagation  (acceleration  data)  using 
segment  model  which  includes  two  throw-off  elements. 


3.3.4  SFEM  Model  for  Force  Identification 

The  SFEM  model  used  in  the  force  identification  process  represents  a  seg¬ 
ment  of  the  beam  structure.  The  implementation  of  the  force  identihcation 
method  presented  here  utilizes  the  response  data  from  the  hve  accelerometer 
positions  in  the  vicinity  of  the  impact  force.  The  spectral  finite  element  model 
is  constructed  by  using  a  throw-off  element  at  each  end  of  the  segment.  For 
these  conditions,  the  frequency-domain  dynamic  stiffness  matrix  is  prepared 
to  relate  the  displacement  and  slope  information  to  the  corresponding  force 
and  moment  information  for  the  beam  segment.  Although  this  model  dif¬ 
fers  considerably  from  the  full  structure  model,  a  response  predicted  by  the 
segment  model  with  the  same  impact  force  previously  considered,  shown  in 
Fig.  4,  matches  perfectly  with  the  local  response  behavior  of  the  full  structure, 
identihed  by  the  box  in  Fig.  2. 

Due  to  the  nature  of  SFEM,  this  force  identihcation  method  calculates 
nodal  forces.  For  conditions  where  the  location  of  the  applied  load  and  the 
accelerometer  are  collocated,  the  force  identihcation  process  does  not  require 
any  additional  steps  in  order  to  obtain  the  identihed  force  information.  How¬ 
ever,  for  cases  where  the  impact  is  applied  between  two  accelerometers,  point 
forces  are  calculated  at  the  neighboring  nodes  when  the  force  identihcation 
procedure  is  applied.  The  summation  of  the  two  force  time  histories  produces 
the  identihed  force  when  certain  conditions  are  satished.  These  conditions  and 
the  justihcation  for  this  process  is  discussed  in  section  3.4  where  the  location 


identification  process  is  presented.  This  approach  allows  for  the  impact  force 
information  to  be  efficiently  identihed  independent  of  the  exact  location  of  the 
impact. 

3.3.5  Numerical  Implementation  and  Parametric  Study 

The  force  identihcation  procedure  is  hrst  applied  to  the  deflection  and  slope 
data  from  the  numerical  simulations.  This  eliminates  the  potential  error  in¬ 
troduced  by  the  deflection  and  slope  calculations  and  allows  for  the  study  to 
focus  on  the  influence  of  other  conditions  on  the  performance  of  the  force  iden¬ 
tihcation  method.  The  theoretical  test  structure  used  to  numerically  study  the 
performance  of  the  force  identihcation  method  is  an  18  ft  (5.49  m)  long  Alu¬ 
minum  beam  with  a  square  cross-section  and  width  of  1  in  (25.4  mm).  Density 
and  Young’s  modulus  values  of  p  =  2700  kg/nY  and  Ey  =  75  GPa  are  used 
and  a  damping  coefficient  value  of  p  =  10^  Ns/m'^  is  identihed  to  match  the 
experimental  response.  The  test  structure  is  chosen  for  geometric  simplicity 
in  order  to  facilitate  experimental  verihcation. 

In  order  to  study  the  ehects  of  various  parameters  on  the  performance  of 
the  force  identihcation  procedure,  an  ideal  condition  is  dehned.  The  theoretical 
test  structure  is  studied  with  free-free  boundary  conditions.  The  full  structure 
simulation  is  conducted  by  using  a  STEM  model  with  66  evenly  spaced  nodes 
and  a  point  impact  force  is  applied.  The  time  window  for  the  numerical 
simulation  is  2.56  s  which  allows  for  the  propagating  mechanical  waves  to 
fully  dissipate  for  the  complete  structure  model  without  throw-oh  elements. 
The  force  information  calculated  with  the  proposed  method  corresponds  to  the 
accelerometer  locations  which  were  included  in  the  subset  of  data  analyzed. 
For  the  ideal  condition,  the  impact  force  is  applied  to  a  node  in  the  middle  of 
the  structure  in  order  to  minimize  the  effect  of  the  reflections  on  the  accuracy 
of  the  identihed  force  information.  The  duration  of  the  impact  is  0.4  ms,  which 
can  be  realized  experimentally. 

The  SEEM  model  and  numerical  simulations  of  the  response  of  the  struc¬ 
ture  are  prepared  and  performed  by  using  MATLAB.  The  proposed  force  iden¬ 
tihcation  method  is  also  implemented  and  studied  by  using  MATLAB  and 
applied  to  the  simulated  response  data  to  study  the  method.  Both  the  simula¬ 
tion  of  the  response  data  and  the  analysis  of  this  data  are  performed  by  using 
SFEM-based  approaches.  This  procedure  is  used  since  other  simulation  tech¬ 
niques,  such  as  time-domain  methods,  have  been  found  to  display  limitations 
with  regard  to  high  frequency  performance  [14].  For  this  reason,  experimental 
response  data  is  also  used  in  order  to  validate  the  performance  of  the  force 
identihcation  method. 
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Figure  5:  Force  information  calculated  for  the  simulated  response  for  an  18  ft 
(5.49  m)  long  beam.  For  this  nominal  case,  there  is  strong  agreement  between 
the  identihed  (connected  points)  and  simulated  (circles)  force  information. 


For  the  ideal  conditions,  the  point  force  is  detected  at  the  location  of  the 
impact,  which  is  at  node  33  at  a  position  of  8.85  ft  (2.70  m)  from  the  reference 
end  of  the  beam.  The  peak  calculated  forces  at  the  other  nodes  within  the 
subset  of  data  analyzed  have  values  less  than  6  N  (<  2%),  which  are  easily 
distinguished  from  the  identihed  impact  force.  The  force  information  identi¬ 
hed  from  the  simulated  response,  shown  in  Fig.  5,  agrees  very  well  with  the 
simulated  force.  In  order  to  quantify  the  accuracy  of  the  force  identihcation 
method,  a  Root  Mean  Square  (RMS)  error  is  calculated  by  comparing  the 
identihed  force  values  to  the  values  used  in  the  simulation.  The  RMS  error 
for  impact  force  identihcation  [27]  provides  quantitative  information  about  the 
average  error  through  the  impact  time  history.  The  quantihcation  of  the  error 
is  performed  by  hrst  separating  the  time  series  into  two  regions:  the  impulsive 
load  and  the  remaining  portion  of  the  time  history,  as  shown  in  Fig.  5.  For  the 
nominal  case,  the  RMS  error  is  2.72  N  during  the  impact  and  an  average  value 
of  2.46  N  exists  for  the  impact  force  considered  over  the  remainder  of  the  1  ms 
presented.  These  errors  are  quite  small  when  compared  with  the  300  N  peak 
value  of  the  impulsive  load,  less  than  1  percent.  The  presence  of  error  in  the 
identihed  force  information  for  the  nominal  case  results  from  the  use  of  the 
segment  model  and  how  it  differs  from  the  complete  structure.  However,  the 
use  of  the  segment  model  provides  the  means  for  the  practical  implementation 
of  the  force  identihcation  method. 


The  geometric  and  loading  properties  corresponding  to  the  results  pre¬ 
sented  in  Fig.  5  provide  nominal  conditions  for  performing  force  characteriza¬ 
tion.  Variations  of  these  properties  are  studied  in  order  to  determine  how  they 
influence  the  ability  of  the  force  identihcation  method  to  accurately  calculate 
the  applied  impulsive  load.  This  includes  the  effects  of  (i)  varying  the  length 
of  the  structure,  (ii)  using  calculated  slope  information  instead  of  simulated 
slope  information,  (iii)  varying  the  location  of  the  applied  force,  (iv)  varying 
the  duration  of  the  impulsive  load,  and  (v)  varying  the  length  of  the  response 
data  analyzed.  During  each  part  of  the  parametric  study,  the  values  of  the  pa¬ 
rameters  not  currently  being  investigated  are  set  to  the  values  corresponding 
to  the  ideal  conditions.  The  acceleration  signals  used  in  the  numerical  studies 
do  not  contain  noise. 
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Figure  6:  Acceleration  response  at  a  node  near  the  center  of  the  beam  for 
(top)  18  ft  (5.49  m)  beam  and  (bottom)  6  ft  (1.83  m)  beam.  Vertical  lines 
identify  when  the  reflection  returns  to  the  impact  location. 


(i)  Structure  Length 

In  order  to  study  the  influence  of  the  length  of  the  structure  on  the  accuracy  of 
the  calculated  force  information,  the  force  identihcation  process  is  applied  to 
response  behavior  simulated  for  a  range  of  lengths  from  0.5  ft  (0.15  m)  to  18  ft 
(5.49  m).  In  each  case,  the  impact  force  is  applied  at  a  node  in  the  middle  of 
the  structure.  While  the  length  of  the  structure  is  varied,  all  other  conditions 
are  maintained  and  the  simulated  displacement  and  slope  information  is  used 
for  the  force  calculation.  The  main  effect  of  varying  the  length  of  the  structure 
is  to  change  the  length  of  time  before  the  rehected  waves  return  the  location 


where  the  impact  force  was  applied.  This  effect  is  illustrated  in  Fig.  6. 

The  vertical  lines  in  Fig.  6  indicate  the  time  when  the  reflection  is  observed 
in  the  acceleration  data.  For  the  nominal  case  with  an  18  ft  (5.49  m)  long  beam, 
the  reflection  is  observed  after  1.71  ms.  When  the  beam  length  is  reduced  to 
6  ft  (1.83  m),  reflections  are  observed  after  only  0.60  ms.  Since  the  shorter 
length  is  one-third  the  initial  length,  it  takes  roughly  1/3  of  the  time  before  the 
reflected  wave  returns.  This  time  depends  on  the  length  of  the  beam  as  well 
as  the  wave  speed.  Since  wave  propagation  is  dispersive  in  beam  structures, 
the  components  of  the  mechanical  waves  at  different  frequencies  travel  with 
different  speeds.  The  return  of  the  reflected  waves  is  identihed  by  detecting 
oscillations  in  the  response  data  after  the  initial  wave. 
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Figure  7:  Force  identihcation  error  with  respect  to  the  length  of  the  structure. 
RMS  error  (*)  during  the  impulsive  load,  average  force  (o)  after  the  impulsive 
load,  and  a  logarithmic  curve  (— )  ht  to  the  RMS  error  values  are  presented. 


Error  information  is  obtained  by  applying  the  force  identihcation  scheme 
to  the  response  data  simulated  for  beam  structures  of  lengths  within  identihed 
range.  The  RMS  error  values  (*)  obtained  by  comparing  the  simulated  and 
identihed  force  values  during  the  impulsive  force  and  the  average  identihed 
force  values  (o)  after  the  impulsive  force  are  presented  in  Fig.  7.  Examination 
of  these  error  values  reveal  that  as  the  length  of  the  structure  is  decreased,  the 
values  of  both  error  measures  appear  to  increase  along  a  logarithmic  prohle. 
This  is  conhrmed  by  using  a  least  squares  approach  to  ht  the  RMS  error  values 
to  a  logarithmic  curve:  Error  =  — 7.55/n  (L)  -|-  14.67  with  E?  =  0.9999.  A 
similar  logarithmic  curve  is  also  found  to  ht  the  trend  of  the  values  of  the  other 
error  measure  (not  shown).  In  addition  to  displaying  the  same  logarithmic 


behavior,  it  can  also  be  seen  in  Fig.  7  that  the  values  of  the  two  error  measures 
are  very  close  at  all  of  the  lengths  considered.  This  suggests  that  the  influence 
of  structure  length  on  the  performance  of  the  force  identihcation  method  is  not 
signihcantly  affected  by  the  magnitude  of  the  force  information.  The  increase 
in  the  error  measure  values  which  are  observed  for  beam  structures  with  shorter 
lengths  is  the  result  of  the  shorter  amount  of  time  before  the  reflected  wave 
returns  to  the  position  on  the  structure  where  the  force  was  applied.  With 
less  time  before  the  reflections  occur,  the  subset  of  data  which  can  be  used  by 
the  force  identihcation  method  is  reduced  and  greater  deviation  exists  between 
the  data  and  the  behavior  expected  by  the  segment  model. 

(ii)  Slope  Calculation 

The  slope  calculation  method  is  tested  by  using  dehection  and  slope  data  from 
numerical  simulations.  The  slope  information  calculated  from  the  simulated 
dehection  data  demonstrates  excellent  agreement  with  the  simulated  slope 
information.  The  maximum  diherence  between  the  calculated  and  simulated 
slope  data  is  between  nine  and  ten  orders  of  magnitude  less  than  the  values  of 
the  slope  from  the  simulation. 

The  error  introduced  by  including  the  slope  calculation  in  the  force  iden¬ 
tihcation  procedure  is  determined  to  be  negligible  for  the  simulated  response 
data.  By  using  the  shape  function  method  to  calculate  the  slope  information, 
the  error  in  the  identihed  impulsive  force  does  not  change  from  the  nominal 
case.  However,  noise  in  the  response  can  cause  the  accuracy  of  the  calculated 
slopes,  and  subsequently  the  identihed  force  information,  to  be  reduced.  The 
calculated  slope  values  and  the  corresponding  identihed  force  information  are 
predicted  to  be  accurate  as  long  as  the  signal-to-noise  ratio  remains  greater 
than  about  100  :  1.  This  limitation  is  addressed  in  the  experimental  verihca- 
tion  by  using  accelerometers  with  sufficiently  high  sensitivity. 

(iii)  Impact  Location 

All  of  the  calculated  force  information  previously  presented  corresponds  to 
impulsive  loading  near  the  center  of  the  structure.  In  order  to  investigate  the 
inhuence  of  the  impact  location  on  the  accuracy  of  the  calculated  force  infor¬ 
mation,  the  impact  force  is  applied  at  a  number  of  locations  on  the  structure. 
The  impact  positions  considered  range  from  1.43  ft  (0.44  m)  to  9.14  ft  (2.79  m) 
from  the  end  of  the  18  ft  (5.49  m)  long  structure.  As  a  result  of  moving  the 
impact  location  closer  to  the  end  of  the  structure,  the  rehections  return  to  the 
impact  location  more  quickly.  The  increase  in  error  resulting  from  moving  the 
impact  location  closer  to  the  end  of  the  structure  is  illustrated  in  Fig.  8. 

The  RMS  error  from  the  identihed  impulsive  load  and  the  average  identihed 
force  values  after  the  impact  both  increase  as  the  impact  position  moves  closer 
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Figure  8:  Force  identification  error  with  respect  to  the  length  of  the  structure. 
RMS  error  (*)  during  the  impulsive  load,  average  force  (o)  after  the  impulsive 
load,  and  a  logarithmic  curve  (— )  fit  to  the  RMS  error  values  are  presented. 


to  the  end  of  the  structure.  Because  moving  the  impact  location  has  the  same 
effect  as  changing  the  length  of  the  structure,  i.e.  reducing  the  length  of  time 
before  the  propagating  wave  returns  to  the  impact  location,  the  error  values 
display  the  same  logarithmic  trend.  In  this  case,  the  logarithmic  curve  fit  to 
the  RMS  error  values  takes  the  form  Error  =  — 4.636L?7,  (ximpact)  +  6.6867 
with  =  0.9916.  When  accounting  for  the  fact  that  when  the  impact  load 
is  applied  in  the  middle  of  the  structure,  the  length  over  which  the  wave 
propagates  before  reflecting  is  half  of  the  length  of  the  structure,  the  data 
in  Fig.  7  and  Fig.  8  are  very  close  to  each  other.  While  all  of  the  conditions 
studied  correspond  to  the  same  geometry  and  material  properties,  these  factors 
influence  the  wave  propagation  speed  and  subsequently  the  length  of  time 
required  for  the  propagating  wave  to  be  reflected  and  return  to  the  impact 
location. 

(iv)  Impact  Duration 

Next,  the  influence  of  impact  duration  on  the  performance  of  the  force  iden¬ 
tification  method  is  studied.  System  responses  are  simulated  for  impulsive 
loads  with  half  and  twice  the  original  duration.  Changing  the  impact  duration 
affects  the  response  of  the  system  in  two  significant  ways.  The  first  and  most 
straightforward  effect  of  changing  the  impact  duration  is  that  the  amount  of 
time  between  the  end  of  the  impulsive  load  and  when  the  reflection  returns 
the  impact  position  is  increased  (reduced  duration)  or  decreased  (extended 
duration).  The  effect  on  the  performance  of  the  force  identification  method 


will  be  similar  to  changing  the  distance  between  the  impact  position  and  the 
end  of  the  structure.  The  second  effect  of  changing  the  impact  duration  is  that 
shorter/longer  duration  loads  will  contain  a  larger /smaller  range  of  frequency 
content.  Since  the  force  identihcation  method  utilizes  the  STEM,  is  it  able  to 
accommodate  the  higher  frequency  content  and  the  change  in  the  error  results 
from  the  travel  time  of  the  propagating  and  then  reflected  wave.  The  two  error 
measures  for  these  conditions  and  the  nominal  case  are  presented  in  Fig.  9. 
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Figure  9:  Force  identihcation  error  with  respect  to  the  impulsive  load  dura¬ 
tion.  RMS  error  (*)  during  the  impulsive  load  and  average  force  (o)  after  the 
impulsive  load  are  presented. 
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From  the  data  presented  in  Fig.  9,  it  can  be  seen  that  for  a  shorter  impact 
duration,  the  RMS  error  remains  the  same  and  the  average  force  after  the 
impact  is  slightly  decreased.  When  the  impact  duration  is  increased  to  0.8  ms, 
both  error  measures  increase  in  value.  The  average  force  after  the  impact 
increases  more  the  RMS  error  during  the  impact.  However,  these  variations 
of  the  impact  duration  result  in  smaller  increases  in  the  error  measures  than 
the  other  parameters. 

(v)  Length  of  Response  Data 

The  influence  of  the  length  of  the  response  data  on  the  accuracy  of  the  re¬ 
constructed  force  is  next  studied.  When  less  of  the  response  data  is  available, 
the  wave  propagation  data  may  not  have  fully  decayed.  This  will  require  the 
application  of  more  aggressive  exponential  windowing  in  order  to  reduce  the 
response  to  zero  by  the  end  of  the  data  set.  When  this  is  done,  the  influence 
of  the  exponential  windowing  on  the  data  at  the  beginning  of  the  data  set 
which  is  used  for  the  force  identihcation  process  is  more  signihcantly  affected. 


This  can  increase  the  error  in  the  identified  force  information.  Properties  such 
as  damping  and  the  boundary  conditions  of  the  structure  can  influence  the 
length  of  time  required  for  the  propagating  wave  to  decay  and  subsequently 
the  length  of  the  data  set  required  to  produce  accurate  results.  In  this  study, 
these  properties  are  kept  constant  and  various  lengths  of  response  data  are 
used  to  calculate  the  applied  force.  The  error  measures  for  the  identified  force 
information  when  using  response  data  sets  of  lengths  ranging  from  5.12  s  to 
down  to  0.64  s  are  presented  in  Fig.  10. 
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Figure  10:  Force  identification  error  with  respect  to  the  length  of  the  response 
data.  RMS  error  (*)  during  the  impulsive  load  and  average  force  (o)  after  the 
impulsive  load  are  presented. 


In  this  figure,  it  can  be  seen  that  above  a  threshold  which  lies  between  0.64  s 
and  1.28  s,  the  error  measures  for  the  identified  force  information  is  close  to 
the  values  for  ideal  conditions.  However,  after  the  length  of  the  response  data 
set  is  decreased  to  0.64  s,  both  of  the  error  measures  increase  by  an  order  of 
magnitude.  These  results  suggest  that  for  a  given  set  of  conditions,  a  certain 
amount  of  time  is  required  for  the  response  of  the  structure  to  sufficiently 
decay  in  order  for  the  force  identification  scheme  to  provide  accurate  results. 
Additionally,  when  the  response  data  of  length  0.64  s  is  used,  the  identified 
force  information  includes  a  consistent  offset  between  the  identified  force  values 
and  the  simulated  force  values,  as  shown  in  Fig.  11.  Based  upon  the  results 
of  this  numerical  study,  the  presence  of  this  characteristic  in  identified  force 
information  indicates  that  a  longer  amount  of  response  data  is  required  to 
improve  the  accuracy  of  the  identified  force. 

Through  this  parametric  study,  the  influence  of  (i)  structure  length,  (ii) 
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Figure  11:  Identified  (connected  points)  and  simulated  (circles)  force  informa¬ 
tion  when  only  one  quarter  of  the  nominal  length  of  the  data  set  is  used. 


slope  calculation,  (iii)  impact  location,  (iv)  impact  duration,  and  (v)  length 
of  response  data  on  the  performance  of  the  force  identihcation  method  was 
studied.  The  accuracy  of  the  reconstructed  impact  force  is  qualitatively  high 
comparing  to  previous  works  [26,  27,  29].  However,  the  variation  in  different 
case  studies  and  measurement  methods  is  a  signihcant  factor  in  the  accuracy 
of  the  results.  For  instance,  using  calculated  slope  information  in  this  work  is 
a  source  of  error  specihc  to  this  work  which  can  be  enhanced  by  using  gyro- 
sensors.  Error  in  the  identihed  force  information  was  determined  to  increase  for 
decreasing  structure  length,  impacts  applied  closer  to  the  end  of  the  structure, 
longer  impact  durations,  and  when  using  response  data  with  lengths  below  a 
system  dependent  threshold.  The  length  of  the  response  data  was  determined 
to  have  the  potential  to  affect  performance  the  most  of  the  parameters  stud¬ 
ied.  Changes  to  the  structure  length,  impact  location,  and  impact  duration 
where  determined  to  all  have  similar  effects  on  the  performance  of  the  force 
identihcation  method.  When  a  signal-to-noise  ratio  of  greater  than  100  :  1 
is  maintained,  the  error  introduced  by  the  slope  calculation  is  expected  to  be 
negligible. 

3.4  Impact  Location  Identification 

The  force  identihcation  method  reconstructs  the  impact  force  in  terms  of  nodal 
forces.  If  the  impact  force  is  applied  between  the  accelerometers,  the  distribu- 


tion  of  the  impact  force  is  utilized  to  locate  the  impact  location.  A  free-body 
diagram  of  the  element  with  the  applied  impact  force  F  and  the  nodal  forces, 
Fi  and  F2,  is  shown  in  Fig.  12.  By  applying  Newton’s  method,  the  summation 
of  vertical  forces  and  the  summation  of  moments  are  calculated  and  presented 
in  Eq.  (38)  and  Eq.  (39),  respectively. 

=  F-Fi-F2  =  mV„  (38) 

=  am  -  Fi—  +  F2—  ^  (^2  ~  ^ 

where  the  mass  of  the  element,  displacement  of  the  center  of  mass  of  the 
element,  and  the  rotation  of  the  element  are  represented  by  m,  Vc,  and  6c, 
respectively.  The  impact  force  is  applied  at  the  distance  x  from  the  hrst  node 
of  the  beam  element  of  length  L. 
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Figure  12:  Free-body  diagram  of  the  beam  element  with  the  impact  force  F 
and  the  nodal  forces,  Fi  and  F2. 


When  the  element  is  relatively  short  and  the  mass  per  length  pA  of  the 
beam  is  small,  the  mass  of  the  element  is  relatively  small  and  the  inertial  force 
can  be  considered  negligible  with  respect  to  the  applied  force.  For  instance, 
the  mass  of  a  3.43  in  (8.71  cm)  element  of  an  Aluminum  beam  with  a  1  in^ 
(6.45  cm^)  cross  section  is  0.137  kg  (weighing  1.34  N)  which  is  quite  small 
compared  to  the  300  N  impact  that  is  applied.  By  considering  the  mass  of 
the  element  negligible,  the  equation  for  the  summation  of  vertical  forces  is 
simplihed  to  Fi  -|-  F2  ~  F. 

When  these  assumptions  are  satished,  the  summation  of  the  two  nodal 
forces  is  equal  to  the  applied  impact  force.  As  a  result,  it  is  possible  to  deter¬ 
mine  the  impact  force  without  knowledge  of  the  exact  location  of  the  impact. 
This  can  signihcantly  reduce  the  computational  costs  when  compared  to  other 
approaches  [27,  28,  29]. 

In  order  to  calculate  the  impact  location  from  the  nodal  forces,  a  straight¬ 
forward  nonlinear  numerical  solver  is  used.  A  flowchart  illustrating  the  algo- 


Figure  13:  Flowchart  of  the  location  identihcation  algorithm. 


rithm  used  for  location  identihcation  is  presented  in  Fig.  13.  The  hrst  step  is 
to  determine  the  distribution  of  the  calculated  impact  force  at  the  two  ends  of 
the  element  by  using  the  experimental  data  collected  from  the  accelerometers. 
The  ratio  of  the  impulses  at  the  two  neighboring  accelerometers  is  used  in 
this  process.  Then,  simulations  are  performed  with  the  beam  model  to  pro¬ 
duce  response  data  for  comparison.  In  order  for  this  process  to  be  effective, 
the  model  must  accurately  represent  the  structure.  Known  material  properties 
and  system  geometry  are  used  for  the  model  and  the  damping  parameter  value 
is  identihed  by  comparing  the  simulated  response  of  the  model  to  the  response 
of  the  experimental  system.  In  the  simulation  portion  of  the  process  illustrated 
in  Fig.  13,  an  extra  node  is  added  between  the  accelerometer  positions  for  ap¬ 
plying  the  impact  force  and  predicting  the  response  of  the  structure.  However, 
only  data  from  the  nodes  which  correspond  to  the  accelerometer  locations  are 
used  in  the  location  identihcation  process.  By  using  the  ratio  of  the  impulses 
for  the  experimental  data  and  the  simulation  data,  the  location  of  the  impact 
is  calculated  with  the  bisection  method. 

Although  the  relationship  between  the  impulse  ratio  and  the  impact  loca¬ 
tion  is  nonlinear,  it  is  also  a  one-to-one  relationship.  This  allows  for  the  use 
of  root-hnding  methods  such  as  the  bisection  method.  In  order  to  start  this 
process,  two  points  are  chosen  close  to  the  two  nodes  of  the  element  on  which 
the  impact  force  was  applied.  The  distance  from  the  two  starting  points  to  the 
nodes  is  chosen  to  be  5%  of  the  length  of  the  element.  This  level  of  precision 
is  selected  as  it  corresponds  to  a  practical  limitation  identihed  in  the  location 
identihcation  procedure.  The  simulated  impulse  ratios  for  the  two  points  are 
then  compared  to  the  impulse  ratio  for  the  experimental  data.  Based  on  the 


Figure  14:  Implementation  of  the  bisection  method  for  the  location  identihca- 
tion  process. 


differences  between  these  ratios,  one  of  the  initial  points  is  discarded  and  a 
simulation  is  performed  to  obtain  a  simulated  impulse  ratio  for  a  new  point 
directly  between  the  two  initial  points.  The  bisection  method  is  continued 
in  this  way,  as  illustrated  in  Fig.  14,  until  the  difference  between  the  experi¬ 
mental  impulse  ratio  and  a  simulated  impulse  ratio  is  less  than  the  prescribed 
tolerance.  The  impact  position  used  for  the  hnal  simulated  case  provides  the 
impact  location  for  the  experimental  data. 

The  ratio  of  the  identihed  impact  forces  at  the  two  nodes  adjacent  to  the 
impact  location,  calculated  by  using  the  simulation  data,  is  shown  in  Fig.  15. 
These  results  are  for  a  specihc  case  and  the  details  of  this  relationship  are 
determined  to  be  influenced  by  material  properties,  element  length,  and  sam¬ 
pling  frequency.  In  order  to  study  the  relationship  between  the  impulse  ratio 
and  the  location  of  the  impact,  the  analytical  model  presented  in  Fig.  12  is 
revisited.  By  using  the  relationship  between  the  impact  force  and  nodal  forces 
and  similarly  considering  the  rotational  inertia  to  be  negligible,  Eq.  (39)  is 
simplihed.  If  the  deformation  of  the  beam  element  is  sufficiently  small  and  the 
moment  difference  AM  is  considered  negligible,  an  analytical  formula  for  the 
impact  force  ratio  is  obtained.  This  formula,  F1/F2  =  {x  —  L)  /x,  provides 
the  approximate  ratio  between  the  two  nodal  forces.  The  impact  force  ratio 
for  these  assumptions  is  plotted  as  a  function  of  position  in  Fig.  15  along  with 
impulse  ratio  data  from  numerical  studies  performed  with  the  SFEM  model  of 
the  beam  structure.  The  analytical  results  qualitatively  agree  with  the  numer¬ 
ical  results.  The  slight  difference  between  the  two  is  believed  to  result  from 
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Figure  15:  Impulse  ratio  as  a  function  of  impact  location  for  the  SFEM  model 
(•)  and  the  analytical  model  (— ).  Red  dashed  lines  indicate  ±  5  %  of  location 
between  nodes 


neglecting  the  deformations  of  the  beam  element  in  the  analytical  model. 

The  impact  force  location  identihcation  process  is  applied  to  simulation 
results  in  order  to  study  the  method.  Representative  results  are  presented 
for  a  case  where  the  impact  force  is  applied  between  nodes  32  and  33  on  the 
18  ft  (5.49  m)  Aluminum  beam  at  20%  of  the  element  length  (0.2L)  from 
node  32.  The  force  information  calculated  from  the  simulated  response  data 
is  presented  in  Fig.  16.  As  anticipated,  larger  force  values  can  be  seen  for 
the  node  closer  to  the  impact  location.  By  combining  the  force  information, 
the  identihed  impact  force  is  calculated.  These  results,  presented  in  Fig.  17, 
display  strong  agreement  between  the  identihed  force  and  the  simulated  force 
without  knowledge  of  the  impact  location.  The  RMS  error  for  the  impact 
and  the  average  for  the  remainder  of  the  2  ms  period  shown  in  the  hgure  are 
2.98  N  and  2.76  N,  respectively.  By  analyzing  the  two  calculated  force  curves 
presented  in  Fig.  16,  an  impulse  ratio  of  8.33  is  calculated.  The  impact  location 
is  calculated  to  be  at  a  distance  of  0.19L  from  node  32.  This  corresponds  to 
an  error  of  1%  of  the  element  length  or  0.034  in  (0.09  cm)  on  the  18  ft  (5.49  m) 
beam  structure. 


Figure  16:  Force  information  calculated  using  simulation  data.  The  force  is 
applied  at  20%  of  the  element  length  from  the  32"''^  node. 


4  Analysis  &  Discussion 

In  this  section,  AWT-FEM  is  used  to  numerically  simulation  nonlinear  wave 
propagation  in  order  to  demonstrate  the  performance  of  the  method.  The 
performance  of  the  force  identihcation  method  is  demonstrated  by  using  ex¬ 
perimentally  collected  data. 

4.1  Modeling  Nonlinear  Wave  Propagation 

In  the  first  subsection,  simulations  on  a  materially  nonlinear  rod  are  conducted 
by  using  AWT-FEM,  AFT-FEM,  and  TFEM.  The  goal  of  this  subsection  is  to 
demonstrate  the  advantages  of  AWT-FEM  in  hdelity  and  computational  effi¬ 
ciency.  Unlike  beams  or  plates,  which  are  strongly  intrinsically  dispersive,  axial 
wave  in  a  rod  has  an  intact  wave  shape  and  it  is  easier  to  identify  the  existence 
of  numerical  errors.  In  the  second  subsection,  simulations  of  a  geometrically 
nonlinear  beam  by  using  AWT-FEM  with  physically  realistic  boundary  con¬ 
ditions  are  conducted.  The  influence  of  the  interaction  with  the  boundaries 
on  the  response  is  analyzed.  In  the  third  subsection,  simulations  of  a  plate 
model  with  a  weak  geometric  nonlinearity  are  presented  to  demonstrate  the 
adaptability  of  the  AWT-FEM  for  2D  models. 


Figure  17:  Summation  of  identified  forces  calculated  using  simulation  data. 
The  force  is  applied  at  20%  of  the  element  length  from  the  32’^'^  node. 


4.1.1  Materially  Nonlinear  Rod 

A  diagram  of  a  hxed-free  rod  subject  to  an  axial  impact  load  at  one  end  is 
shown  in  Fig.  18.  The  parameters  of  the  system  are  presented  in  Table.  2. 
The  rod  is  modeled  with  a  nonlinear  constitutive  relationship  as  defined  in 
Eqn.  (16)  with  a  coefficient  a  =  20. 

Simulation  results  for  the  materially  nonlinear  rod  model  in  Eqn.  (17)  with 
AWT-FEM,  AFT-FEM,  and  TFEM  are  shown  in  Fig.  19.  The  Daubechies 
wavelet  with  an  order  of  iV  =  14  is  used  for  the  spectrally-uncoupled  wavelet 
transform  and  can  yield  sufficient  smoothness  in  the  responses.  Fifty  elements 
are  used  for  the  AWT-FEM  and  the  AFT-FEM  and  350  elements  are  used 
for  the  TFEM.  For  systems  with  physically  realistic  boundaries,  additional 
damping  combined  with  a  longer  time-window  is  required  for  the  AFT-FEM 
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Figure  19:  Velocity  responses  at  the  impacted  position  on  the  rod  obtained  by 
using  AWT-FEM,  AFT-FEM,  and  TFEM. 


to  avoid  wrap-around  errors.  However,  convergence  issues  can  occur  for  the 
AFT-FEM  iterative  procedure  with  long  time-window.  In  the  simulation,  the 
AFT-FEM  is  not  able  to  produce  a  reliable  physical  response.  In  order  to 
obtain  a  converged  result  for  comparison  with  other  methods,  a  large  damping 
factor  T]  =  1  X  combined  with  the  original  time  window  is  used. 

An  impulse  in  the  velocity  is  created  by  the  impact  load  at  0.1  ms.  The 
wave  is  reflected  by  the  hxed  end  and  travels  back  to  the  impacted  position 
at  0.9  ms.  The  wave  shape  in  the  reflected  wave  is  shifted  to  the  left-hand 
side.  This  distortion  of  wave  shape  caused  by  the  nonlinearity  is  capture 
by  all  these  three  methods.  Due  to  the  additional  damping,  the  amplitude 
response  simulated  by  using  AFT-FEM  is  severely  reduced.  Wrap-around 
errors  also  occurs  in  the  time-history.  These  results  demonstrate  that  the  AFT- 
FEM  cannot  work  for  nonlinear  structures  without  the  use  of  non-reflecting 
boundaries.  In  order  to  leak  the  energy  out  of  the  system,  either  a  long 
time  window  combined  with  artihcial  damping  or  a  non-reflecting  boundary 
at  one  end  is  needed.  The  high  hdelity  performance  of  the  AFT-FEM  is  only 
possible  with  the  help  of  these  two  mechanisms.  In  the  comparison,  a  physical 
boundary  is  presented  and  the  same  time-window  is  adopted  for  all  three 
numerical  methods.  The  performance  of  AFT-FEM  is  greatly  compromised. 
Due  to  the  periodic  nature  of  FFT,  the  traveling  signal  gets  wrapped  back 
into  the  reflected  wave  and  results  in  periodic  errors  in  the  response.  The 
TFEM  can  produce  a  response  with  adequate  resolution  at  the  expense  of  a 
rehned  meshing  (350  elements).  However,  in  the  insert  of  Fig.  19,  the  reflected 


Table  2:  System  parameters  of  rod. 
Parameter  Value 

70  GPa 


Elastic  modulus,  E 
Cross  section,  A 
Mass  density,  p 
Rod  length,  L 
Time  window,  T 
Impact  duration,  Tp 
Impact  amplitude,  Ta 
Sampling  frequency,  / 
Nonlinear  coefficient,  a 


TT  X  25  mmx25  mm 
2800  kg/m^ 

2  m 
1  ms 
50  ps 
100  kN 
1000  kHz 
20 


Table  3:  Computation  time. 


AWT-FEM 

AFT-FEM 

TFEM 

11.9  s 

6.3  s 

25.8  s 

wave  is  affected  by  spurious  oscillations  late  in  the  time  series  resulting  from 
erroneously  introduced  numerical  dispersion.  For  an  impact  load  with  shorter 
duration  and  higher  frequency  content,  this  error  can  be  further  magnihed 
and  completely  distort  the  response  [14].  On  the  contrary,  the  AWT-FEM 
only  uses  50  elements  to  produce  an  intact  wave  response  and  capture  the 
nonlinear  distortion  with  high  hdelity.  The  computation  time  for  these  three 
methods  is  listed  in  Table.  3.  A  Dell  desktop  with  a  Intel  Core  Quad  CPU  is 
used  for  all  simulations.  The  parallel  computing  toolbox  in  MAT  LAB  is  used 
to  execute  the  calculation  at  different  wavelet  point  in  parallel.  With  the  use  of 
parallel  computing,  the  AWT-FEM  only  uses  11.9  s  with  50  elements  to  obtain 
an  accurate  response  while  the  TFEM  needs  25.8  s  with  350  elements.  The 
AFT-FEM  cannot  produce  a  correct  response  for  these  conditions.  However, 
its  fast  computational  performance  indicates  that  if  a  fast  numerical  technique 
similar  to  FFT  can  be  applied  to  modify  the  wavelet  transform  procedure,  the 
computational  performance  of  the  AWT-FEM  will  be  further  improved. 

A  convergence  study  of  the  AWT-FEM  and  the  TFEM  is  conducted.  The 
same  impact  condition  in  Table.  2  is  adopted.  Simulations  start  with  1  ele¬ 
ments  and  incrementally  rehne  the  mesh  by  adding  10  elements.  The  error 
is  dehned  as  the  absolute  value  of  the  relative  difference  between  the  current 
simulation  result  of  the  velocity  response  Vcur  and  the  reference  in  the  pre¬ 
vious  state  Vpre  at  the  impacted  position,  as  dehned  in  Eqn.  (40.  The  state 
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Figure  20:  Comparison  of  convergence  between  the  AWT-FEM  and  the  non¬ 
linear  TFEM  when  the  impact  amplitude  Ta  =  100  kN  and  the  nonlinear 
coefficient  a  =  20. 


being  considered  can  be  number  of  elements  or  order  of  Daubechies  wavelet 
depending  on  the  convergence  study.) 

n 

E\\VcuriU)-Vpre{U)V\\ 

Error  =  — - - - .  (40) 

E  \\Vcur{U)V\\ 
i=l 

As  show  in  Fig.  20,  the  convergence  performance  of  the  AWT-FEM  is 
signihcantly  better  than  the  TFEM.  With  50  elements,  the  AWT-FEM  has  a 
4.2%  error  level  while  the  error  of  TFEM  is  68.8%.  By  increasing  the  number 
of  elements  to  90,  the  AWT-FEM  can  reach  a  0.9%  error  level  while  the  TFEM 
still  has  an  error  of  23.1%. 

A  convergence  study  of  the  AWT-FEM  with  respect  to  the  order  of  Daubechies 
wavelet  is  also  conducted  with  the  same  error  measure  dehned  in  Eqn.  (40). 

As  shown  in  Fig.  21,  the  order  of  Daubechies  wavelet  starts  from  4  to  20  with 
an  increment  of  2.  When  the  order  is  greater  and  equal  to  8,  the  error  level  is 
below  1%. 

4.2  Geometrically  Nonlinear  Beam 

A  diagram  of  a  cantilevered  beam  with  clamped-free  boundaries  is  studied  as 
shown  in  Fig.  22.  A  point  impact  load  with  amplitude  150  kN  is  applied  at  the 
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Figure  21:  Convergence  of  the  order  of  Daubechies  wavelet. 


free  end.  The  parameters  of  the  system  are  presented  in  Table  4.  Simulations 
for  a  geometrically  nonlinear  beam  model  in  Eqn.  (21)  and  Eqn.  (22)  with 
the  clamped-free  boundaries  is  conducted  by  using  AWT-FEM.  A  Daubechies 
wavelet  with  an  order  of  A^  =  14  is  used  for  the  spectrally-uncoupled  wavelet 
transform  and  can  yield  sufficient  smoothness  in  the  response.  Fifty  elements 
are  used  in  the  spectral  element  formulation  to  provide  satisfactory  accuracy 
and  sufficient  resolution. 

A  comparison  of  the  velocity  response  at  the  impacted  position  obtained  us¬ 
ing  the  AWT-FEM  and  a  nonlinear  TFEM  is  shown  in  Fig.  23.  The  AWT-FEM 
uses  50  elements  and  the  TFEM  uses  100  elements.  As  reported  in  previous  re¬ 
search,  the  geometric  nonlinearities  have  limited  influence  on  transverse  wave 
propagation  in  beams  in  the  intermediate  strain  regime  [71,  72]  compared  to 
the  nonlinear  dispersion  in  the  rod  model.  Since  the  beam  is  a  strong  intrinsi¬ 
cally  nonlinear  system,  the  influence  of  the  nonlinear  dispersion  introduced  by 


t  (ms) 


Figure  23:  Velocity  responses  at  the  impacted  position  on  the  beam  obtained 
by  using  AWT-FEM  and  TFEM. 


geometric  nonlinearity  on  the  intrinsic  dispersion  relationship  of  the  beam  in 
the  intermediate  strain  regime  is  negligible.  The  nonlinear  responses  obtained 
using  two  methods  are  close  to  each  other.  The  difference  in  the  later  part  of 
the  response,  where  the  reflected  signal  interacts  with  the  forward  traveling 
waves,  may  be  attributed  to  the  extra  dispersion  introduced  by  the  TFEM. 

A  convergence  study  of  the  AWT-FEM  and  the  TFEM  for  the  beam  model 
is  conducted.  The  same  impact  condition  in  Table.  4  is  adopted.  Simulations 
start  with  1  and  10  elements  and  incrementally  rehne  the  mesh  by  adding  10 
elements.  The  error  measure  dehned  in  Eqn.  (40)  is  adopted  here.  The  error 
level  of  the  AWT-FEM  is  two  orders  of  magnitude  smaller  than  the  TFEM. 
With  50  elements,  the  AWT-FEM  has  only  0.1%  error  while  the  TFEM  has 
3.8%.  By  increasing  the  number  of  elements  to  90,  the  error  of  the  AWT-FEM 
reduces  to  0.04%  while  the  TFEM  is  0.2%. 

A  convergence  study  of  the  AWT-FEM  with  respect  to  the  order  of  Daubechies 
wavelet  is  also  conducted  with  the  same  error  measure  dehned  in  Eqn.  (40)  for 
the  beam  model.  As  shown  in  Fig.  25,  the  order  of  Daubechies  wavelet  starts 
from  4  to  20  with  an  increment  of  2.  When  the  order  is  greater  than  and  equal 
to  6,  the  error  level  is  below  0.5%. 

4.2.1  Geometrically  Nonlinear  Plate 

A  diagram  of  a  cantilevered  plate  with  clamped-free  boundaries  along  the  x 
dimension  and  free-free  boundaries  along  the  y  dimension  is  shown  in  Fig.  26. 

A  point  impact  load  is  applied  at  the  middle  point  of  edge  BC  {Ly).  Other 


Table  4:  System  parameters  of  a  beam 


Parameter 

Value 

Elastic  modulus,  E 

70  GPa 

Cross  section,  A 

25  mm  x25  mm 

Moment  of  inertia,  I 

3.26  X  10-®  m^ 

Mass  density,  p 

2800  kg/m^ 

Beam  length,  L 

1  m 

Time  window,  T 

1  ms 

Impact  duration,  Tp 

50  /is 

Impact  amplitude,  Ta 

100  kN 

Sampling  frequency,  / 

1000  kHz 

Table  5:  System  parameters  of  a  plate 


Parameter 

Value 

Elastic  modulus,  E 

70  GPa 

Density,  p 

2800  kg/m^ 

Length  1, 

1  m 

Length  2,  Ly 

0.1  m 

Thickness,  h 

25  mm 

Poisson’s  ratio,  u 

0.33 

Time  window,  T 

1  ms 

Impact  amplitude,  Em 

100  kN 

Impact  duration,  Tp 

100  ps 

Sampling  frequency,  / 

500  kHz 

general  distributed  load  can  also  be  chosen.  The  parameters  of  the  system 
are  presented  in  Table.  5.  Along  the  y  dimension,  free-free  boundaries  have  no 
constraints  on  the  axial  motion.  Along  the  x  dimension,  the  length  is  chosen  to 
be  La;  =  1  m,  the  same  as  was  used  for  the  beam  model  in  previous  subsection. 
The  width  is  chosen  to  be  L^,  =  0.1  m  to  approximate  a  narrow  2D  cantilever 
beam. 

Simulations  for  a  plate  model  with  a  weak  geometric  nonlinearity  based  on 
Eqn.  (31)  in  the  intermediate  strain  value  regime  is  conducted  by  using  AWT- 
FEM.  A  Daubechies  wavelet  with  an  order  of  A^  =  22  is  used  for  the  wavelet 
transform  with  respect  to  time.  A  Daubechies  wavelet  with  an  order  of  A^  =  16 
with  a  sampling  rate  Ay  =  0.04  m  is  used  for  the  wavelet  transform  with 
respect  to  the  spatial  coordinate  y.  These  wavelets  are  able  to  provide  sufficient 
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Figure  24:  Comparison  of  convergence  between  the  AWT-FEM  and  the  non¬ 
linear  TFEM  when  the  impact  amplitude  Ta  =  100  kN 


smoothness  in  the  response.  Twenty  elements  are  used  to  approximate  both 
the  X  dimension  and  the  y  dimension. 

The  comparison  of  the  velocity  responses  at  the  impacted  position  obtained 
using  the  AWT-FEM,  ID  TFEM,  and  2D  TFEM  are  shown  in  Fig.  27.  The 
ID  TFEM  with  100  elements  is  applied  onto  a  ID  beam  model  with  the  same 
length  as  in  the  previous  subsection.  Since  the  width  of  the  plate  is  narrow, 
the  response  of  the  AWT-FEM  has  a  similar  trend  with  the  response  of  the  ID 
TFEM.  There  are  three  major  differences  between  them.  First,  the  width  of 
the  plate  makes  it  stiffer  than  the  corresponding  ID  beam  model,  which  results 
in  a  smaller  amplitude  of  the  response  in  the  AWT-FEM.  Second,  the  response 
in  the  AWT-FEM  after  the  initial  pulse  includes  small  amplitude  oscillations 
due  to  the  interaction  with  reflected  waves  from  both  sides  while  the  response 
in  the  ID  TFEM  is  smooth.  Third,  the  later  part  of  the  responses  in  the 
two  methods  when  the  reflected  waves  travel  back  is  different.  This  may  be 
related  to  the  effect  of  the  sides  of  the  plate  and  the  influence  of  the  different 
stiffnesses  on  the  wave  speed.  A  2D  TFEM  with  a  50  x  20  elements  using 
AN SYS  is  also  adopted  here  to  produce  the  wave  propagation  process.  As 
shown  by  the  dashed-doted  green  line,  the  response  is  signihcantly  different 
from  the  other  two.  For  this  problem  with  extreme  impact  loading,  the  2D 
TFEM  with  a  reasonable  hne  mesh  using  a  sequential  time  integrator  cannot 
produce  a  convergent  result.  It  is  also  worth  noting  that  by  further  increasing 
the  impact  load  and  leading  the  strain  into  a  strong  value  regime,  all  three 


Figure  25:  Convergence  of  the  order  of  Daubechies  wavelet. 


methods  will  encounter  convergence  issue. 

4.3  Force  Identification 

The  performance  of  the  force  identihcation  method  is  further  studied  with 
experimental  response  data  in  order  to  verify  the  performance  observed  by 
using  simulated  response  data.  The  experimental  setup  is  a  6  ft  (1.83  m)  long 
aluminum  bar  with  a  square  cross-section  and  width  of  1  in  (25.4  mm),  as 
shown  in  Fig.  28.  It  is  shorter  than  the  ideal  beam  structure  dehned  in  the 
numerical  study  in  order  to  accommodate  laboratory  space  limitations.  The 
structure  is  instrumented  with  twenty-two  PCB  356A22  accelerometers,  with 
a  frequency  range  of  0.5  to  4000  Hz,  distributed  evenly  along  its  length.  This 
number  of  accelerometers  allows  for  good  resolution  in  the  slope  information 
calculated  from  the  acceleration  data.  These  are  triaxial  accelerometers  (axial, 
transverse  and  lateral  directions),  but  only  transverse  response  data  is  used. 


t  (ms) 


Figure  27:  Velocity  responses  at  the  impacted  position  on  the  beam  obtained 
using  AWT-FEM,  ID  TFEM,  and  2D  TFEM. 


The  acceleration  data  is  collected  with  a  sampling  frequency  of  102.4  kHz 
in  order  to  capture  the  response  with  good  temporal  resolution.  The  data  is 
down-sampled  by  a  factor  of  four  in  order  to  accommodate  memory  limitations 
when  applying  the  force  identihcation  procedure.  The  structure  is  suspended 
from  bungee  cords  in  order  to  approximate  free-free  boundary  conditions.  A 
PCB  086D05  modal  impact  hammer  is  used  to  apply  impulsive  loads  in  the 
vertical  direction.  The  applied  force  is  measured  by  the  force  transducer  in 
the  impact  hammer  for  comparison  with  the  identihed  force  information.  For 
transducer  signal  conditioning  and  A/D  signal  conversion,  a  LMS  SCADAS 
III  acquisition  system  is  used.  The  measured  data  is  recorded  by  using  LMS 
Test. Lab  Software,  Ver.  IIA.  The  acceleration  data  sets  are  chosen  from  the 
accelerometers  in  the  vicinity  of  the  impact  location.  The  acceleration  data 
contains  the  response  signal  of  the  structure  due  to  the  impact  and  also  the 
reflections  due  to  the  boundaries  of  the  structure.  The  acceleration  data  is 
dispersive  due  to  the  dispersive  nature  of  the  structure. 

The  impact  force  identihcation  method  is  applied  to  experimental  response 
data  collected  for  conditions  when  the  impact  location  is  collocated  with  an 
accelerometer  and  when  the  impact  is  applied  between  accelerometers.  Rep¬ 
resentative  results  for  these  conditions  are  presented  and  discussed  below. 

As  it  was  shown  in  the  parametric  study,  the  calculated  impact  force  is 
more  precise  at  cases  where  the  impact  occurs  close  to  the  middle  of  the 
structure.  The  hrst  data  set  presented  corresponds  to  an  impact  force  being 
applied  to  the  experimental  system  at  the  location  of  accelerometer  number  11 


Figure  28:  Photograph  of  the  experimental  setup. 


near  the  middle  of  the  structure.  The  force  information  obtained  by  applying 
the  force  identihcation  method  is  presented  in  Fig.  29.  Since  the  impact  force 
was  applied  at  the  location  of  the  accelerometer,  the  location  of  the  impact 
is  easily  identihed  as  the  position  of  accelerometer  11.  The  identihed  force 
information  agrees  well  with  the  measured  force,  successfully  capturing  the 
qualitative  characteristics  of  the  impulsive  load.  However,  some  discrepancies 
are  observed.  This  error  is  believed  to  be  introduced  as  a  result  of  errors  in 
the  slope  calculation  due  to  signal  noise  and  also  due  to  the  low  frequency 
stop-band  of  the  accelerometers  which  are  used  to  perform  modal  analysis. 
The  RMS  error  for  this  case  is  34.60  N  and  an  average  value  of  10.45  N 
exists  during  the  remainder  of  the  1  ms  time  window  presented  in  Fig.  29. 
In  comparison  with  the  ideal  conditions,  the  length  of  the  beam  is  shorter 
and  noise  is  present  in  the  acceleration  data.  While  each  of  these  differences 
contributes  to  the  reduced  accuracy  of  the  calculated  force  information,  the 
RMS  error  for  the  impulsive  load  is  only  12%  of  the  maximum  force  value. 

The  next  two  data  sets  presented  were  collected  when  the  impact  force 
was  applied  between  accelerometers.  These  data  sets  are  used  in  order  to  test 
the  force  identihcation  as  well  as  the  location  identihcation  methods.  The 
impact  force  identihcation  method  is  applied  to  calculate  the  forces  applied 
to  the  structure  and  the  location  identihcation  method  is  used  to  determine 
where  the  impact  force  was  applied.  The  hrst  of  these  two  data  sets  corre¬ 
spond  to  conditions  where  the  impact  force  is  applied  at  the  center  between 
accelerometers  11  and  12.  The  impact  force  information  obtained  from  the 
force  identihcation  method  is  presented  in  Fig.  30.  The  detected  force  values 
at  nodes  11  and  12  exhibit  similar  properties  since  the  impact  was  applied 
at  equal  distances  from  the  two  accelerometers.  However,  after  the  impulsive 
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Figure  29:  Force  information  calculated  from  the  experimental  data,  force 
applied  at  node  11. 


load,  the  identihed  force  values  increase  or  decrease  linearly  with  time. 

The  summation  of  the  identihed  forces  at  the  accelerometers  on  either 
side  of  the  impact  location  provides  the  identihed  force.  The  identihed  force 
information  (A)  agrees  well  with  the  measured  force  information  (o)  during 
the  impact.  After  the  impact,  the  calculated  force  exhibits  a  persistent  ohset 
from  the  measured  force  values.  The  RMS  error  for  the  summation  of  the  two 
calculated  forces  is  39.20  N  during  the  impact  and  averages  53.21  N  during 
the  remainder  of  the  1  ms.  Based  on  the  parametric  study,  the  force  ohset 
after  the  impact  is  believed  to  be  due  to  the  short  length  of  the  acceleration 
data  set.  Access  to  a  longer  set  of  the  response  data  is  expected  to  yield  more 
accurate  results.  The  bungee  cords  supporting  the  experimental  system  also 
apply  a  small  but  constant  force  to  the  structure  which  is  not  addressed  in  the 
model  and  may  inhuence  the  results.  By  using  the  impulse  ratio  calculated  for 
the  force  values  at  accelerometers  11  and  12,  the  impact  location  is  calculated 
to  be  at  a  distance  of  0.52L  from  accelerometer  11.  This  corresponds  to  a 
deviation  of  only  2%  of  the  distance  between  the  accelerometers  which  is  less 
than  2  mm. 

The  hnal  data  set  correspond  to  conditions  where  the  impact  force  is  ap¬ 
plied  at  a  distance  of  0.40L  from  accelerometer  11.  The  impact  force  informa¬ 
tion  obtained  from  the  force  identihcation  method  is  presented  in  Fig.  31.  As 
the  numerical  study  predicted,  the  calculated  force  values  for  the  accelerom¬ 
eter  nearest  to  the  impact  location  are  larger  than  the  values  for  the  other 
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Figure  30:  Force  information  calculated  using  experimental  data  for  the  force 
applied  at  50%  of  the  distance  between  the  accelerometers  from  accelerometer 
11. 


location.  Again,  the  calculated  force  values  at  the  positions  of  the  two  ac¬ 
celerometers  gradually  increase  or  decrease  after  the  impact.  The  identihed 
force  information  is  obtained  from  the  summation  of  the  calculated  forces  and 
it  agrees  well  with  the  measured  force  during  the  impact.  The  RMS  error  is 
35.61  N  during  the  impact  and  averages  83.78  N  during  the  remainder  of  the 
1  ms  presented  in  Fig.  31.  The  impact  force  is  identihed  to  be  at  a  distance 
of  0.37L  from  accelerometer  11.  This  corresponds  to  a  deviation  of  only  3% 
of  the  distance  between  the  accelerometers  which  is  less  than  3  mm  from  the 
actual  location  of  the  applied  impact  force.  Considering  the  fact  that  the  im¬ 
pact  was  manually  applied  and  the  hammer  tip  is  6  mm  wide,  the  errors  in 
the  results  of  the  location  identihcation  method  are  acceptable. 

The  experimental  verihcation  illustrates  a  good  agreement  between  the 
experiment  and  the  simulation  predictions.  However,  the  force  identihcation 
method  of  the  experimental  setup  might  experience  errors  due  to  the  discrep¬ 
ancies  in  the  material  properties,  loading  conditions,  mass  of  accelerometers 
and  other  diherent  unknown  factors.  These  errors  can  be  reduced  by  conduct¬ 
ing  diherent  tests  in  order  to  calibrate  the  method. 
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Figure  31:  Force  information  calculated  using  experimental  data  for  the  force 
applied  at  40%  of  the  distance  between  the  accelerometers  from  accelerometer 
11. 

5  Closure 

In  this  work,  numerical  simulation  and  force  identihcation  &  localization  meth¬ 
ods  were  developed  to  facilitate  the  high-hdelity  study  of  nonlinear  wave  propa¬ 
gation  caused  by  extreme  impact  loading  conditions.  By  using  a  wavelet-based 
spectral  hnite  element  method  and  an  alternating  wavelet-time  framework,  the 
linear  method  was  adapted  for  nonlinear  systems.  The  alternating  wavelet¬ 
time  hnite  element  method  (AWT-FEM)  is  demonstrated  to  provide  high- 
hdelity  results,  be  compatible  with  parallel  computing,  and  provide  a  means 
to  study  general  nonlinear  behavior  in  a  range  of  simple  structures.  The  direct 
iteration  algorithm  adopted  in  this  study  limits  the  convergence  performance 
of  the  AWT-FEM.  For  ID  beam  and  2D  plate  with  geometric  nonlinearity, 
the  AWT-FEM  is  able  to  obtain  a  converged  response  within  an  intermediate 
strain  value  regime  (10“®  to  10“^)  where  the  inhuence  of  the  geometric  nonlin¬ 
earity  is  limited.  In  order  to  address  this  issue,  a  modihed  Newton-Raphson 
iterative  method  was  subsequently  developed  to  replace  the  direct  iteration 
algorithm  in  order  to  improve  the  convergence  performance  of  the  AWT-FEM 
and  facilitate  its  application  for  large  deformation  problems. 

An  impact  force  identihcation  method  using  the  spectral  hnite  element 
method  was  also  presented  and  demonstrated  with  a  beam  structure.  With  this 
method,  the  impact  force  can  be  determined  without  precise  information  about 
the  location  of  the  impact  on  the  structure.  The  procedure  was  demonstrated 


with  simulated  response  data  for  propagating  mechanical  waves  and  validated 
by  using  experimental  data.  In  simulations,  excellent  agreement  was  observed 
for  nominal  conditions.  Sources  of  error  in  the  identihed  force  information  were 
investigated  through  a  parametric  study.  The  most  influential  parameter  in 
the  force  calculation  procedure  is  the  length  of  the  response  data  set.  Reduced 
performance  was  identihed  when  decreasing  the  structure  length,  moving  the 
impact  location  toward  the  end  of  the  structure,  and  increasing  the  impact 
duration.  Also,  the  trends  associated  with  varying  the  structure  length,  impact 
position,  and  impact  duration  were  studied.  The  force  identihcation  method 
was  also  applied  to  experimental  data  and  was  able  to  successfully  identify  the 
characteristics  of  the  impact  and  provide  identihed  impact  force  information 
with  relatively  low  error.  In  conjunction  with  the  force  identihcation  method, 
a  new  technique  for  accurately  determining  the  location  of  the  impact  was 
proposed.  When  applied  to  experimental  data,  the  impact  locations  were 
identihed  within  5%  of  the  distance  between  the  accelerometers  from  the  actual 
location,  which  corresponds  to  an  error  of  less  than  0.17  in  (0.4  cm)  on  a  6  ft 
(1.83  m)  structure. 

The  AWT-FEM  is  promising  tool  with  great  potential  for  studying  a  wide 
range  of  nonlinear  wave  propagation  in  various  structural  components.  While 
the  numerical  studies  conducted  in  this  work  provided  great  insight  into  the 
performance  of  the  new  numerical  simulation  technique,  future  studies  will 
incorporate  signihcant  experimental  verihcation  in  order  to  further  study  and 
rehne  the  tool.  Future  work  on  the  new  methods  for  force  identihcation  and 
localization  will  explore  the  ohset  in  the  calculate  force  information  which  fol¬ 
lows  the  impulsive  load.  Further  rehnements  to  the  force  identihcation  meth¬ 
ods  will  also  include  incorporating  more  advanced  sensor  technology,  such  as 
gyroscopic  sensors  to  directly  measure  rotational  response  information.  While 
these  force  identihcation  methods  have  been  demonstrated  for  beam  struc¬ 
tures,  its  application  to  structures  such  as  rods,  plates,  and  more  complicated 
structures  such  as  shells,  trusses  and  arches  will  also  be  explored. 
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